Project

General

Profile

makeEffHist_v3.C

Version 3 - Muhammad Elnimr, 02/17/2015 12:56 PM

 
1
#define makeEffHist_cxx
2
#include "makeEffHist.h"
3
#include <TH2.h>
4
#include <TStyle.h>
5
#include <TCanvas.h>
6

    
7
void makeEffHist::Loop(int nEvents=0,char * filePrefex="")
8
{
9
//   In a ROOT session, you can do:
10
//      Root > .L makeEffHist.C
11
//      Root > makeEffHist t
12
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
13
//      Root > t.Show();       // Show values of entry 12
14
//      Root > t.Show(16);     // Read and show values of entry 16
15
//      Root > t.Loop();       // Loop on all entries
16
//
17
  prettyPlots();
18
//     This is the loop skeleton where:
19
//    jentry is the global entry number in the chain
20
//    ientry is the entry number in the current Tree
21
//  Note that the argument to GetEntry must be:
22
//    jentry for TChain::GetEntry
23
//    ientry for TTree::GetEntry and TBranch::GetEntry
24
//
25
//       To read only selected branches, Insert statements like:
26
// METHOD1:
27
//    fChain->SetBranchStatus("*",0);  // disable all branches
28
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
29
// METHOD2: replace line
30
//    fChain->GetEntry(jentry);       //read all branches
31
//by  b_branchname->GetEntry(ientry); //read only this branch
32
   if (fChain == 0) return;
33

    
34
  double fMinMCKE=0.05;           // Minimum MC particle kinetic energy (GeV).
35
  double fMinMCLen=20.;          // Minimum MC particle length in tpc (cm).
36
  double fMatchColinearity=0.97;  // Minimum matching colinearity.
37
  double fMatchDisp=2.0;         // Maximum matching displacement.
38
  double fWMatchDisp=15;        // Maximum matching displacement in the w direction.
39
  double pdg=-13;
40

    
41
   Long64_t nentries = fChain->GetEntriesFast();
42
   TH1F * ntracks=new TH1F("ntracks","No. of Tracks",100,0,10);
43
   TH2F * trkytrkz=new TH2F("trkytrkz","Track Spts. y vs. z",1000,-10,300,1000,-150,150);
44
   TH2F * trkytrkx=new TH2F("trkytrkx","Track Spts. y vs. x",1000,-50,250,1000,-150,150);
45
   TH2F * trkztrkx=new TH2F("trkztrkx","Track Spts. z vs. x",1000,-50,250,1000,-100,200);
46
   TH3F * trkztrkxtrky=new TH3F("xyz","xyz",10,-10,300,10,-100,200,10,-10,300);
47

    
48

    
49
   TH2F * hgoodlen=new TH2F("hgoodlen","hgoodlen",300,0,300,300,0,300);
50

    
51
   TH1F * htrkx=new TH1F("htrkx","Track Spacepoints x (cm)",100,-100,200);
52
   TH1F * htrky=new TH1F("htrky","Track Spacepoints y (cm)",100,-100,200);
53
   TH1F * htrkz=new TH1F("htrkz","Track Spacepoints z (cm)",100,-100,200);
54
   TH1F * htrkstartx=new TH1F("htrkstartx","Track Start x (cm)",100,-100,300);
55
   TH1F * htrkstarty=new TH1F("htrkstarty","Track Start y (cm)",100,-100,200);
56
   TH1F * htrkstartz=new TH1F("htrkstartz","Track Start z (cm)",100,-100,200);
57
   TH1F * htrkstartx_diff=new TH1F("htrkstartx_diff","htrkstartx_diff",100,-50,50);
58
   TH1F * htrkstarty_diff=new TH1F("htrkstarty_diff","htrkstarty_diff",100,-50,50);
59
   TH1F * htrkstartz_diff=new TH1F("htrkstartz_diff","htrkstartz_diff",100,-50,50);
60
   TH1F * htrkendx=new TH1F("htrkendx","Track End x (cm)",100,-100,300);
61
   TH1F * htrkendy=new TH1F("htrkendy","Track End y (cm)",100,-100,200);
62
   TH1F * htrkendz=new TH1F("htrkendz","Track End z (cm)",100,-100,200);
63
   TH1F * hntracks_reco=new TH1F("hntracks_reco","No. reco tracks",10,-1,4);
64

    
65
   TH1F * htrkspacepoints=new TH1F("htrkspacepoints","All spacepoints",50,0,100);
66
   TH1F * hgtrkspacepoints=new TH1F("hgtrkspacepoints","Good spacepoints",50,0,100);
67
   TH1F * hetrkspacepoints=new TH1F("hetrkspacepoints","Eff. spacepoints",50,0,100);
68

    
69
   TH1F * htrkclust=new TH1F("htrkclust","All clust",50,0,100);
70
   TH1F * hgtrkclust=new TH1F("hgtrkclust","Good clust",50,0,100);
71
   TH1F * hetrkclust=new TH1F("hetrkclust","Eff. clust",50,0,100);
72

    
73
   TH1F * htrkhits=new TH1F("htrkhits","All hits",25,0,400);
74
   TH1F * hgtrkhits=new TH1F("hgtrkhits","Good hits",25,0,400);
75
   TH1F * hetrkhits=new TH1F("hetrkhits","Eff. hits",25,0,400);
76

    
77

    
78
   TH2F * hnHitsnClust=new TH2F("hnHitsnClust","nhits vs. nclusts",300,-10,290,50,-2,48);
79
   TH2F * hnHitsnSpcpts=new TH2F("hnHitsnSpcpts","nhits vs. nspcpts",300,-10,290,100,-10,90);
80
   TH2F * hnClustSpcpts=new TH2F("hnClustSpcpts","nclusts vs. nspcpts",50,-2,48,100,-10,90);
81

    
82
   Long64_t nbytes = 0, nb = 0;
83
   MCHists mchists("somename");
84
   if(nEvents==0);
85
   else
86
     nentries=nEvents;
87
   int COUNTER=0;
88
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
89
      Long64_t ientry = LoadTree(jentry);
90
      if (ientry < 0) break;
91
      nb = fChain->GetEntry(jentry);   nbytes += nb;
92
      // if (Cut(ientry) < 0) continue;
93
      double pdg=-13;
94
      int  nTracks=0;
95
      int nspacepoints=0;
96
      int skipEvent=0;
97
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
98
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
99
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
100
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101
      //
102
      //
103
      //  Reco track loop
104
      //
105
      //
106
      int totalSpacepoints=0;
107
      for(int j=0;j<ntracks_reco;j++)
108
        {
109
          /*          if(trkpdg[j]!=-13 )
110
            {
111
              cout << "----" << trkpdg[j] << endl;
112
              cout << " track stuff " << trkd2[j] << endl;
113
              }*/
114
          //          if(trklen_cut_MC[0]<5  && trklen[j]>60 ) std::cout << " Event No." <<  event << std::endl;
115
          if(!(1/*trkang[j]<0 &&trkd2[j]<0.01 &&trkmatchdisp[j]>0  && abs(trklen[j]-10-trklen_cut_MC[0]+10)<30 */))
116
            {
117
            skipEvent=1;
118
            break;
119
            }
120
          totalSpacepoints+=ntrkhits[j];
121
        }
122
      /*      if(trklen_cut_MC[0]<200)
123
              continue;*/
124
      if(trklen_cut_MC[0]==0)
125
        continue;
126
      if(skipEvent) 
127
        continue;
128
      /*      if(ntracks_reco==0)
129
              continue;  */
130
      if(trkpdg_MC[0]!=-13)
131
        continue;
132
      
133
      //      if(trklen[j]-10-trklen_cut_MC[0]+10<30)
134
      //      if(trklen_cut_MC[0]>180 &&trklen[i]< 100 )
135
      //        cout << "Event number is :" << event << endl;
136

    
137
      /*      mchists.fHmcstartx->Fill(trkstartx_MC[0]);
138
      mchists.fHmcstarty->Fill(trkstarty_MC[0]);
139
      mchists.fHmcstartz->Fill(trkstartz_MC[0]);
140
      mchists.fHmcendx->Fill(trkendx_MC[0]);
141
      mchists.fHmcendy->Fill(trkendy_MC[0]);
142
      mchists.fHmcendz->Fill(trkendz_MC[0]);*/
143
      mchists.fHmctheta->Fill(trktheta_MC[0]);
144
      mchists.fHmcphi->Fill(trkphi_MC[0]);
145
      mchists.fHmctheta_xz->Fill(trktheta_xz_MC[0]);
146
      mchists.fHmctheta_yz->Fill(trktheta_yz_MC[0]);
147
      mchists.fHmcmom->Fill(trkmom_MC[0]);
148
      mchists.fHmclen->Fill(trklen_cut_MC[0]);
149
      mchists.fHmclenshort->Fill(trklen_cut_MC[0]);
150
      hntracks_reco->Fill(ntracks_reco);
151
      htrkspacepoints->Fill(totalSpacepoints);
152
      htrkclust->Fill(nclust);
153

    
154
      int nhits_trk=0;
155
      /*      for(int p=0;p<nhits;p++)
156
        {
157
          if(hit_trkid[p]>0)
158
            nhits_trk++;
159
            }*/
160
      //      htrkhits->Fill(nhits_trk);
161
      htrkhits->Fill(nhits);
162

    
163
      COUNTER++;
164
      
165
      
166
      if(ntracks_reco==0 )
167
        {
168
          if(trklen_cut_MC[0]>200)
169
            {
170
              cout << "=============================================" << endl; 
171
              cout <<" Entry : " << jentry << " Run : " << run << "Subrun :  "  << subrun << " Event : "  << event <<endl;
172
              cout << "Zero Tracks Reconstructed" << endl;
173
              cout << "=============================================" << endl; 
174
            }
175
        }
176
      else
177
        {
178
          bool foundmatch = false;
179
          for(int j=0;j<ntracks_reco;j++)
180
            {
181
              if( trklen[j]/trklen_cut_MC[0] > 0.5)
182
                {
183
                  //                  if(trklen_cut_MC[0]>200)
184
                    {
185
                      foundmatch = true;
186
                    }
187
                }
188
            }
189
          if (!foundmatch){
190
            cout << "=============================================" << endl;
191
            cout << " Run : " << run << "Subrun :  "  << subrun << " Event : "  << event <<endl;
192
            cout << "Bad event : MC length (within readoutwindow) " << trklen_cut_MC[0]<<endl;
193
            cout << "=============================================" << endl; 
194
          }
195
          /*
196
          for(int j=0;j<ntracks_reco;j++)
197
            {
198
              if(!( trklen[j]/trklen_cut_MC[0] > 0.5))
199
                {
200
                  if(trklen_cut_MC[0]>200)
201
                    {
202
                      cout << "=============================================" << endl; 
203
                      cout << " Entry : " << jentry << " Run : " << run << "Subrun :  "  << subrun << " Event : "  << event <<endl;
204
                      cout << "Bad event : Reco length " << trklen[j] <<" MC length (within readoutwindow) " << trklen_cut_MC[0]<<endl;
205
                      cout << "=============================================" << endl; 
206
                    }
207
                }
208
                }
209
          */
210
        }
211
      
212
      
213

    
214
      for(int j=0;j<ntracks_reco;j++)
215
        {
216
          /*          if( 0 
217
             //trkd2[j]>0.2
218
               || trklen_cut_MC[0]<fMinMCLen
219
               || trkcolinearity[j]<fMatchColinearity 
220
              // ||trklen[j]/trklen_cut_MC[0] < 0.5
221
               || abs(trkwmatchdisp[j]) <= fWMatchDisp &&  trkmatchdisp[j] <= fMatchDisp 
222
               ||ntrkhits[j]==0
223
               )
224
               continue;*/
225
          //          if(abs(trklen[j]-trklen_cut_MC[j]+100)>10 ) continue;
226
          //if(abs(trklen[j]-trklen_cut_MC[j])>10 ) continue;
227
          //          if(ntrkhits[j]==0) continue;
228
          //          cout << "Event no. is : " << event << endl;
229
          /*          htrkstartx->Fill(trkstartx[j]);
230
          htrkstarty->Fill(trkstarty[j]);
231
          htrkstartz->Fill(trkstartz[j]);*/
232

    
233
        
234
          htrkendx->Fill(trkendx[j]);
235
          htrkendy->Fill(trkendy[j]);
236
          htrkendz->Fill(trkendz[j]);
237

    
238
          
239

    
240
          htrkstartx->Fill(trkstartx[j]);
241
          htrkstarty->Fill(trkstarty[j]);
242
          htrkstartz->Fill(trkstartz[j]);
243

    
244
          htrkstartx_diff->Fill(trkendx[j]-trkstartx[j]);
245
          htrkstarty_diff->Fill(trkendy[j]-trkstarty[j]);
246
          htrkstartz_diff->Fill(trkendz[j]-trkstartz[j]);
247

    
248
          mchists.fHstartdx->Fill(trkstartx_MC[0]-trkstartx[j]);
249
          mchists.fHstartdy->Fill(trkstarty_MC[0]-trkstarty[j]);
250
          mchists.fHstartdz->Fill(trkstartz_MC[0]-trkstartz[j]);
251
          for(int tt=0;tt<ntrkhits[j];tt++)
252
            {
253
              
254
              trkytrkz->Fill(trkz[j][tt],trky[j][tt]);
255
              trkytrkx->Fill(trkx[j][tt],trky[j][tt]);
256
              trkztrkx->Fill(trkx[j][tt],trkz[j][tt]);
257
              trkztrkxtrky->Fill(trkz[j][tt],trky[j][tt],trkx[j][tt]);
258
              
259
              htrkx->Fill(trkx[j][tt]);
260
              htrky->Fill(trky[j][tt]);
261
              htrkz->Fill(trkz[j][tt]);
262
              
263
            }
264
          
265
          nspacepoints+=ntrkhits[j];
266
          nTracks++;
267
          mchists.fHtrkd2->Fill(trkd2[j]);
268
          //          if( trkmatchdisp[j] <= fMatchDisp  && trklen_cut_MC[0]> -999)
269
          //          if(  trklen_cut_MC[0]> -999)
270
          {
271
            mchists.fHlvsl->Fill(trklen_cut_MC[0],trklen[j]);
272
            mchists.fHlvslvstrkang->Fill(trkang[j],trklen_cut_MC[0],trklen[j]);
273
            mchists.fHdl->Fill(trklen[j]-trklen_cut_MC[0]);
274
          }
275

    
276
          bool good = (/*trkang[j]>0 &&ntracks_reco==1&& trklen_cut_MC[0]>10 &&abs(trkwmatchdisp[j]) <= fWMatchDisp &&  trkmatchdisp[j] <= fMatchDisp &&*/
277
                       //                        trklenratio[j] > 0.5 
278
                       trklen[j]/trklen_cut_MC[0] > 0.5
279
                       //                                && pdg==-13 
280
                        //                        && trklen_cut_MC[0]>fMinMCLen
281
                        //&& trkcolinearity[j]>fMatchColinearity
282
                        ) ;
283
          //          bool good=1;
284
          //trklen[j]/trklen_cut_MC[0] > 0.5 && pdg==13 && trklen_cut_MC[0]>fMinMCLen  && trkcolinearity[j]>fMatchColinearity);
285
          if(good) {
286
            /* mchists.fHgstartx->Fill(mcstart.X());
287
               mchists.fHgstarty->Fill(mcstart.Y());
288
            mchists.fHgstartz->Fill(mcstart.Z());
289
            mchists.fHgendx->Fill(mcend.X());
290
            mchists.fHgendy->Fill(mcend.Y());
291
            mchists.fHgendz->Fill(mcend.Z());*/
292
            mchists.fHgtheta->Fill(trktheta_MC[j]);
293
            mchists.fHgphi->Fill(trkphi_MC[j]);
294
            mchists.fHgtheta_xz->Fill(trktheta_xz_MC[j]);
295
            mchists.fHgtheta_yz->Fill(trktheta_yz_MC[j]);
296
            mchists.fHgmom->Fill(trkmom_MC[j]);
297
            hgoodlen->Fill(trklen_cut_MC[0],trklen[j]);
298
            hgtrkspacepoints->Fill(totalSpacepoints);
299
            mchists.fHglen->Fill(trklen_cut_MC[j]);
300
            mchists.fHglenshort->Fill(trklen_cut_MC[j]);
301
            //            hgtrkhits->Fill(nhits_trk);
302
            hgtrkhits->Fill(nhits);
303
            hgtrkclust->Fill(nclust);
304
          }
305
          else{
306
           
307
          }
308
        }
309

    
310
      //      hnHitsnClust->Fill((double)nhits_trk,(double)nclust);
311
      //      hnHitsnSpcpts->Fill((double)nhits_trk,(double)totalSpacepoints);
312
      hnHitsnClust->Fill((double)nhits,(double)nclust);
313
      hnHitsnSpcpts->Fill((double)nhits,(double)totalSpacepoints);
314
      hnClustSpcpts->Fill((double)totalSpacepoints,(double)nclust);
315
      //
316
      //
317
      // End of  Reco track loop
318
      //
319
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
320
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
321
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
322
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
323
      //
324
      ntracks->Fill(nTracks);
325
      //      cout << " event no " << jentry<< " no of spacepoints " << nspacepoints << endl;
326
   }
327

    
328
   cout <<"+++++++++++ "<< endl;
329
   cout << "::effcalc" << endl;
330
   cout <<"+++++++++++ "<< endl;
331
   effcalc(mchists.fHgtheta, mchists.fHmctheta, mchists.fHetheta);
332
   cout <<"+++++++++++ "<< endl;
333
   effcalc(mchists.fHgphi, mchists.fHmcphi, mchists.fHephi);
334
   cout <<"+++++++++++ "<< endl;
335
   effcalc(mchists.fHgtheta_xz, mchists.fHmctheta_xz, mchists.fHetheta_xz);
336
   cout <<"+++++++++++ "<< endl;
337
   effcalc(mchists.fHgtheta_yz, mchists.fHmctheta_yz, mchists.fHetheta_yz);
338
   cout <<"+++++++++++ "<< endl;
339
   effcalc(mchists.fHgmom, mchists.fHmcmom, mchists.fHemom);
340
   cout <<"+++++++++++ "<< endl;
341
   effcalc(mchists.fHglen, mchists.fHmclen, mchists.fHelen);
342
   cout <<"+++++++++++ "<< endl;
343
   effcalc(mchists.fHglenshort, mchists.fHmclenshort, mchists.fHelenshort);
344
   cout <<"+++++++++++ "<< endl;
345
   effcalc( hgtrkspacepoints,  htrkspacepoints, hetrkspacepoints);
346
   cout <<"+++++++++++ "<< endl;
347
   effcalc( hgtrkhits,  htrkhits, hetrkhits);
348
   cout <<"+++++++++++ "<< endl;
349
   effcalc( hgtrkclust,  htrkclust, hetrkclust);
350

    
351

    
352
   //   cout << COUNTER << endl;
353
   //
354
   //
355
   // theta
356
   //
357
     
358
   /*   new TCanvas;
359
   mchists.fHmctheta->Draw();
360
   new TCanvas;
361
   mchists.fHgtheta->Draw();*/
362
  
363
   cout<< "+++++++++++" << endl;
364
   cout << " Phi " << endl;
365
   cout<< "+++++++++++" << endl;
366
    new TCanvas;
367
   mchists.fHmcphi->Draw();
368
   mchists.fHmcphi->SetFillColor(kPink);
369
   //   new TCanvas;
370
   mchists.fHgphi->Draw("same");
371
   mchists.fHgphi->SetFillColor(kBlue);
372
   cout<< "+++++++++++" << endl;
373
   cout << " theta " << endl;
374
   cout<< "+++++++++++" << endl;
375
   new TCanvas;
376
   mchists.fHmctheta->Draw();
377
   mchists.fHmctheta->SetFillColor(kPink);
378
   //   new TCanvas;
379
   mchists.fHgtheta->Draw("same");
380
   mchists.fHgtheta->SetFillColor(kBlue);
381
   cout<< "+++++++++++" << endl;
382
   cout << " Momentum " << endl;
383
   cout<< "+++++++++++" << endl;
384
   new TCanvas;
385
   mchists.fHmcmom->Draw();
386
   mchists.fHmcmom->SetFillColor(kPink);
387
   //  new TCanvas;
388
   mchists.fHgmom->Draw("same");
389
   mchists.fHgmom->SetFillColor(kBlue);
390
   cout<< "+++++++++++" << endl;
391
   cout << " Length " << endl;
392
   cout<< "+++++++++++" << endl;
393
   new TCanvas;
394
   mchists.fHmclen->Draw();
395
   mchists.fHmclen->SetMaximum(700);
396
   mchists.fHmclen->SetMinimum(0);
397
   mchists.fHmclen->SetFillColor(kPink);
398
   //  new TCanvas;
399
   mchists.fHglen->Draw("same");
400
   mchists.fHglen->SetFillColor(kBlue);
401
   //   mchists.fHglen->SetFillStyle(3445);
402
   /////////
403
   /////////
404
   /////////
405
   /////////
406
   /////////
407

    
408
  TString tmpName="display"; 
409
  tmpName+="_hist.root";
410
  TFile * file=new TFile(tmpName,"UPDATE");
411

    
412
   TString tmpSt="";
413
   new TCanvas;
414
   mchists.fHephi->Draw("p");
415
   tmpSt=  mchists.fHephi->GetName();
416
   tmpSt+=filePrefex;
417
   mchists.fHephi->SetName(tmpSt);
418
   mchists.fHephi->Write();
419
   new TCanvas;
420
   mchists.fHetheta->Draw("p");
421
   tmpSt=  mchists.fHetheta->GetName();
422
   tmpSt+=filePrefex;
423
   mchists.fHetheta->SetName(tmpSt);
424
   mchists.fHetheta->Write();
425
   new TCanvas;
426
   mchists.fHetheta_xz->Draw("p");
427
    tmpSt=  mchists.fHetheta_xz->GetName();
428
   tmpSt+=filePrefex;
429
   mchists.fHetheta_xz->SetName(tmpSt);
430
   mchists.fHetheta_xz->Write();
431
   new TCanvas;
432
   mchists.fHetheta_yz->Draw("p");
433
  tmpSt=  mchists.fHetheta_yz->GetName();
434
   tmpSt+=filePrefex;
435
   mchists.fHetheta_yz->SetName(tmpSt);
436
   mchists.fHetheta_yz->Write();
437
   new TCanvas;
438
   mchists.fHelen->Draw("p");
439
     tmpSt=  mchists.fHelen->GetName();
440
   tmpSt+=filePrefex;
441
   mchists.fHelen->SetName(tmpSt);
442
   mchists.fHelen->Write();
443
   new TCanvas;
444
   mchists.fHelenshort->Draw("p");
445
   tmpSt=  mchists.fHelenshort->GetName();
446
   tmpSt+=filePrefex;
447
   mchists.fHelenshort->SetName(tmpSt);
448
   mchists.fHelenshort->Write();
449

    
450
   new TCanvas;
451
   mchists.fHlvsl->Draw();
452
  tmpSt=  mchists.fHlvsl->GetName();
453
   tmpSt+=filePrefex;
454
   mchists.fHlvsl->SetName(tmpSt);
455
   mchists.fHlvsl->Write();
456
   //c1->SaveAs("lvsl.png");
457
   new TCanvas;
458
   mchists.fHdl->Draw();
459
  tmpSt=  mchists.fHdl->GetName();
460
   tmpSt+=filePrefex;
461
   mchists.fHdl->SetName(tmpSt);
462
   mchists.fHdl->Write();
463
   //c1_n2->SaveAs("dl.png");
464
   new TCanvas;
465
   mchists.fHlvslvstrkang->Draw();
466
   new TCanvas;
467
   mchists.fHemom->Draw();
468
   tmpSt=  mchists.fHemom->GetName();
469
   tmpSt+=filePrefex;
470
   mchists.fHemom->SetName(tmpSt);
471
   mchists.fHemom->Write();
472

    
473
   
474
   new TCanvas;
475
   mchists.fHmctheta_xz->Draw();
476
   mchists.fHgtheta_xz->Draw("same");
477
   mchists.fHgtheta_xz->SetFillColor(kGray);
478
   mchists.fHgtheta_xz->SetMinimum(0);
479
   new TCanvas;
480
   mchists.fHmctheta_yz->Draw();
481
   mchists.fHgtheta_yz->Draw("same");
482
   mchists.fHgtheta_yz->SetFillColor(kGray);
483
   
484
   new TCanvas;
485
   hgoodlen->Draw("zcol");
486
   
487
   new TCanvas;
488
  htrkstartx_diff->Draw();
489
  htrkstartx_diff->SetFillColor(kGray);
490
  
491
  new TCanvas;
492
  htrkstarty_diff->Draw();
493
  htrkstarty_diff->SetFillColor(kGray);
494
  
495
  new TCanvas;
496
  htrkstartz_diff->Draw();
497
  htrkstartz_diff->SetFillColor(kGray);
498
  
499
  new TCanvas;
500
  hetrkspacepoints->Draw("p");
501
  tmpSt=  hetrkspacepoints->GetName();
502
  tmpSt+=filePrefex;
503
  hetrkspacepoints->SetName(tmpSt);
504
  hetrkspacepoints->Write();
505

    
506
  new TCanvas;
507
  hetrkhits->Draw("p");
508
  tmpSt=  hetrkhits->GetName();
509
  tmpSt+=filePrefex;
510
  hetrkhits->SetName(tmpSt);
511
  hetrkhits->Write();
512

    
513
 new TCanvas;
514
  hetrkclust->Draw("p");
515
  tmpSt=  hetrkclust->GetName();
516
  tmpSt+=filePrefex;
517
  hetrkclust->SetName(tmpSt);
518
  hetrkclust->Write();
519

    
520
  new TCanvas;
521
  hnHitsnClust->Draw("zcol");
522
  hnHitsnClust->Write();
523
  new TCanvas;
524
  hnHitsnSpcpts->Draw("zcol");
525
  hnHitsnSpcpts->Write();
526
  new TCanvas;
527
  hnClustSpcpts->Draw("zcol");
528
  hnClustSpcpts->Write();
529
  file->Close();
530
  //  hgtrkspacepoints->Draw("same");
531
  //htrkspacepoints->Draw();
532
  
533
  //  new TCanvas;
534
  //  mchists.htrkd2->Draw();
535
  //  new TCanvas
536
  //   saveHist(mchists.fHtrkd2);
537
  // new TCanvas;
538
  //   mchists.fHstartdz->Draw();
539
  // saveHist( mchists.fHdl, mchists.fHlvsl);
540
  //  saveHist(filePrefex,trkytrkz,trkytrkx,trkztrkx,htrkstartx,htrkstarty,htrkstartz,mchists.fHdl,mchists.fHlvsl,htrkendx,htrkendy,htrkendz,hntracks_reco);
541
  /*            mchists.fHmcmom->Draw();*/
542

    
543
}
544

    
545
void effcalc(const TH1* hnum, const TH1* hden, TH1* heff)
546
{
547
  int nbins = hnum->GetNbinsX();
548
  if (nbins != hden->GetNbinsX())
549
    cout << "something is wrong" << endl;
550
  if (nbins != heff->GetNbinsX())
551
    cout << "something is wrong" << endl;
552
  // Loop over bins, including underflow and overflow.
553
  cout << " Hist Name:    "<< heff->GetName()  << endl;
554
  for(int ibin = 0; ibin <= nbins+1; ++ibin) {
555
    double num = hnum->GetBinContent(ibin);
556
    double den = hden->GetBinContent(ibin);
557
    if(den == 0.) {
558
      heff->SetBinContent(ibin, 0.);
559
      heff->SetBinError(ibin, 0.);
560
    }
561
    else {
562
      double eff = num / den;
563
      if(eff < 0.)
564
        eff = 0.;
565
      /*      if(eff > 1.)
566
              eff = 1.;*/
567
      double err = std::sqrt(eff * (1.-eff) / den);
568
      //double err = std::sqrt(num*num+den*den);
569
      heff->SetBinContent(ibin, eff);
570
      heff->SetBinError(ibin, err);
571
    }
572
    if( heff->GetBinContent(ibin)>0.8)
573
      {
574
        cout << " Bin No" << ibin << endl;
575

    
576
        cout << "+++ Before" <<  endl;
577
        cout << " Bin Content" << heff->GetBinContent(ibin-3) << endl;
578
        cout << " Bin Error" << heff->GetBinError(ibin-3) << endl;
579
        cout << " num" << hnum->GetBinContent(ibin-3) << endl;
580
        cout << " den" << hden->GetBinContent(ibin-3) << endl;
581

    
582
        cout << " Bin Content" << heff->GetBinContent(ibin-2) << endl;
583
        cout << " Bin Error" << heff->GetBinError(ibin-2) << endl;
584
        cout << " num" << hnum->GetBinContent(ibin-2) << endl;
585
        cout << " den" << hden->GetBinContent(ibin-2) << endl;
586

    
587
        cout << " Bin Content" << heff->GetBinContent(ibin-1) << endl;
588
        cout << " Bin Error" << heff->GetBinError(ibin-1) << endl;
589
        cout << " num" << hnum->GetBinContent(ibin-1) << endl;
590
        cout << " den" << hden->GetBinContent(ibin-1) << endl;
591

    
592
        cout << "++++++" << endl;
593
        cout << " Bin Content" << heff->GetBinContent(ibin) << endl;
594
        cout << " Bin Error" << heff->GetBinError(ibin) << endl;
595
        cout << " num" << hnum->GetBinContent(ibin) << endl;
596
        cout << " den" << hden->GetBinContent(ibin) << endl;
597
      }
598

    
599
  }
600
  heff->SetMinimum(0.);
601
  heff->SetMaximum(1.05);
602
  heff->SetMarkerStyle(20);
603
}
604

    
605
void saveHist(char * filePrefex,TH2F * h1,TH2F *h2,TH2F* h3,TH1F * h4, TH1F * h5, TH1F * h6,TH1F * h7, TH2F *h8,TH1F * h9, TH1F * h10, TH1F * h11, TH1F * h12)
606
{
607
 
608

    
609

    
610
  TCanvas * tmpCan=new TCanvas("tmpCan","tmpCan",1200,800);
611
  TString tmpSt="";
612
  tmpCan->Divide(4,3);
613
  tmpCan->cd(1);
614
  h1->Draw();
615
  tmpSt=h1->GetName();
616
  tmpSt+=filePrefex;
617
  h1->SetName(tmpSt);
618
  h1->Write();
619
  h1->SetFillColor(kGray);
620
  tmpCan->cd(2);
621
  h2->Draw();
622
  tmpSt=h2->GetName();
623
  tmpSt+=filePrefex;
624
  h2->SetName(tmpSt);
625
  h2->SetFillColor(kGray);
626
  h2->Write();
627
  tmpCan->cd(3);
628
  h3->Draw();
629
 tmpSt=h3->GetName();
630
  tmpSt+=filePrefex;
631
  h3->SetName(tmpSt);
632
  h3->SetFillColor(kGray);
633
  h3->Write();
634
  tmpCan->cd(4);
635
  h4->Draw();
636
 tmpSt=h4->GetName();
637
  tmpSt+=filePrefex;
638
  h4->SetName(tmpSt);
639
  h4->SetFillColor(kGray);
640
  h4->Write();
641
  tmpCan->cd(5);
642
  h5->Draw();
643
  h5->SetFillColor(kGray);
644
  h5->Write();
645
  tmpCan->cd(6);
646
  h6->Draw();
647
 tmpSt=h6->GetName();
648
  tmpSt+=filePrefex;
649
  h6->SetName(tmpSt);
650
  h6->SetFillColor(kGray);
651
  h6->Write();
652
  tmpCan->cd(7);
653
  h7->Draw();
654
 tmpSt=h7->GetName();
655
  tmpSt+=filePrefex;
656
  h7->SetName(tmpSt);
657
  h7->SetFillColor(kGray);
658
  h7->Write();
659
  tmpCan->cd(8);
660
  h8->Draw("zcol");
661
 tmpSt=h8->GetName();
662
  tmpSt+=filePrefex;
663
  h8->SetName(tmpSt);
664
  h8->SetFillColor(kGray);
665
  h8->Write();
666
  tmpCan->cd(9);
667
  h9->Draw();
668
 tmpSt=h9->GetName();
669
  tmpSt+=filePrefex;
670
  h9->SetName(tmpSt);
671
  h9->SetFillColor(kGray);
672
  h9->Write();
673
  tmpCan->cd(10);
674
  h10->Draw();
675
   tmpSt=h10->GetName();
676
  tmpSt+=filePrefex;
677
  h10->SetName(tmpSt);
678
  h10->SetFillColor(kGray);
679
  h10->Write();
680
  tmpCan->cd(11);
681
  h11->Draw();
682
   tmpSt=h11->GetName();
683
  tmpSt+=filePrefex;
684
  h11->SetName(tmpSt);
685
  h11->SetFillColor(kGray);
686
  h11->Write();
687
  tmpCan->cd(12);
688
  h12->Draw();
689
 tmpSt=h12->GetName();
690
  tmpSt+=filePrefex;
691
  h12->SetName(tmpSt);
692
  h12->SetFillColor(kGray);
693
  h12->Write();
694
  tmpCan->Draw();
695

    
696

    
697
  //  tmpCan->SaveAs(histName);
698
  //  delete tmpCan;
699
  //  tmpCan->SaveAs(histName);
700
  //  cout << histName << endl;
701
}
702
void saveHist(TH1F * h1, TH2F * h2)
703
{
704
  TCanvas * tmpCan=new TCanvas("tmpCan","tmpCan",600,800);
705
  tmpCan->Divide(1,2);
706
  tmpCan->cd(1);
707
  h1->Draw();
708
  tmpCan->cd(2);
709
  h2->Draw();
710
  tmpCan->Draw();
711
  //  tmpCan->SaveAs(histName);
712
  //  delete tmpCan;
713
  //  tmpCan->SaveAs(histName);
714
  //  cout << histName << endl;
715
}
716

    
717

    
718
void prettyPlots()
719
{
720
  //    gROOT  -> SetStyle("Default");
721
  //    gROOT  -> SetStyle("Plain");
722
  gStyle->SetNdivisions(110,"xyz");
723
  gStyle -> SetStatFont(42);
724
  //  gStyle -> SetOptStat(111);
725
  gStyle -> SetOptStat(111);
726
  gStyle -> SetOptFit(1111);
727
  gROOT  ->ForceStyle();
728
  gStyle -> SetPalette(1);
729
  gStyle->SetStatX(0.8201149);
730
  gStyle->SetStatY(0.8444915);
731
  //  gStyle->SetStatW(0.8807471);
732
  //  gStyle->SetStatH(0.8644068);
733

    
734
  gStyle -> SetLabelFont(42,"XYZ");
735
  gStyle -> SetTitleFont(42, "XYZ");
736
  gStyle -> SetLabelSize(0.04, "XYZ");
737
  gStyle -> SetTitleSize(0.05, "XYZ");
738
  gStyle -> SetTitleOffset(1.0, "X");
739
  gStyle -> SetTitleOffset(1.0, "Y");
740
  gStyle -> SetTitleOffset(1.0, "Z");
741
}