]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TRD/Macros/AliTRDseedsAnalysis.C
A few mods in the GRP preprocessor (to get some greens finally)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDseedsAnalysis.C
1 #ifndef __CINT__
2   #include <iostream.h>
3   #include "AliTRDtracker.h"
4   #include "AliTRDcluster.h" 
5   #include "AliTRDv1.h"
6   #include "AliTRDgeometry.h"    
7   #include "AliTRDparameter.h"    
8   #include "alles.h"  
9   #include "AliTRDmcTrack.h"
10   #include "AliTRDtrack.h"
11   
12   #include "TFile.h"
13   #include "TParticle.h"
14   #include "TStopwatch.h"
15 #endif
16
17
18 void AliTRDseedsAnalysis() {
19
20   const Int_t nPrimaries = (Int_t) 3*86030/4;
21   //  const Int_t nPrimaries = 500;
22
23   Bool_t page[6]; for(Int_t i=0; i<6; i++) page[i]=kTRUE;
24
25   const Int_t nPtSteps = 30;
26   const Float_t maxPt = 3.;
27   const Float_t minPt = 0.;
28   
29   const Int_t nEtaSteps = 40;
30   const Float_t maxEta = 1.;
31   const Float_t minEta = -1.;
32   
33   // page[0]
34   TH1F *hNcl=new TH1F("hNcl","Seeds No of Clusters",255,-9.5,500.5); 
35   hNcl->SetFillColor(4);
36   hNcl->SetXTitle("No of clusters"); 
37   hNcl->SetYTitle("counts"); 
38   TH1F *hNtb=new TH1F("hNtb","Seeds No of TB with clusters",160,-9.5,150.5); 
39   hNtb->SetFillColor(4);
40   hNtb->SetXTitle("No of timebins with clusters"); 
41   hNtb->SetYTitle("counts"); 
42   TH1F *hNep=new TH1F("hNep","Seeds end point",160, -9.5, 150.5); 
43   hNep->SetFillColor(4);
44   hNep->SetXTitle("outermost timebin with cluster"); 
45   hNep->SetYTitle("counts"); 
46   TH1F *hNmg=new TH1F("hNmg","Seeds max gap",160, -9.5, 150.5); 
47   hNmg->SetFillColor(4);
48   hNmg->SetXTitle("max gap (tb)"); 
49   hNmg->SetYTitle("counts"); 
50
51   // page[1]
52   Float_t fMin = -5.5, fMax = 134.5;
53   Int_t   iChan = 140; 
54   TH2F *h2ep = new TH2F("h2ep","MC vs RT (end point)",iChan,fMin,fMax,iChan,fMin,fMax); 
55   h2ep->SetMarkerColor(4);
56   h2ep->SetMarkerSize(2);
57   TH2F *h2ntb = new TH2F("h2ntb","MC vs RT (TB with Clusters)",iChan,fMin,fMax,iChan,fMin,fMax); 
58   h2ntb->SetMarkerColor(4);
59   h2ntb->SetMarkerSize(2);
60   TH2F *h2mg = new TH2F("h2mg","MC vs RT (max gap)",iChan/2,fMin,fMax/2,iChan/2,fMin,fMax/2); 
61   h2mg->SetMarkerColor(4);
62   h2mg->SetMarkerSize(2);
63
64   // page[2]
65   TH1F *hPt_all=new TH1F("hPt_all","Seeds Pt",nPtSteps,minPt,maxPt);
66   hPt_all->SetLineColor(4);
67   hPt_all->SetXTitle("Pt (GeV/c)"); 
68   hPt_all->SetYTitle("counts"); 
69
70   TH1F *hPt_short=new TH1F("hPt_short","Short seeds Pt",nPtSteps,minPt,maxPt);
71   hPt_short->SetLineColor(8);
72   TH1F *hPt_long=new TH1F("hPt_long","Long seeds Pt",nPtSteps,minPt,maxPt);
73   hPt_long->SetLineColor(2);
74   
75   TH1F *hrtPt_short=new TH1F("hrtPt_short","RT from short seeds",nPtSteps,minPt,maxPt);
76   hrtPt_short->SetLineColor(8);
77   TH1F *hrtPt_long=new TH1F("hrtPt_long","RT from long seeds",nPtSteps,minPt,maxPt);
78   hrtPt_long->SetLineColor(2);
79
80   TH1F *hEta_all=new TH1F("hEta_all","Seeds Eta",nEtaSteps,minEta,maxEta);
81   hEta_all->SetLineColor(4);
82   hEta_all->SetXTitle("Eta"); 
83   hEta_all->SetYTitle("counts"); 
84   TH1F *hEta_short=new TH1F("hEta_short","Short seeds Eta",nEtaSteps,minEta,maxEta);
85   hEta_short->SetLineColor(8);
86   TH1F *hEta_long=new TH1F("hEta_long","Long seeds Eta",nEtaSteps,minEta,maxEta);
87   hEta_long->SetLineColor(2);
88   
89   TH1F *hrtEta_short=new TH1F("hrtEta_short","RT from short seeds",nEtaSteps,minEta,maxEta);
90   hrtEta_short->SetLineColor(8);
91   TH1F *hrtEta_long=new TH1F("hrtEta_long","RT from long seeds",nEtaSteps,minEta,maxEta);
92   hrtPt_long->SetLineColor(2);
93
94   // page[3]
95   TH1F *hEff_long = new TH1F("hEff_long","Efficiency vs Pt",nPtSteps,minPt,maxPt); 
96   hEff_long->SetLineColor(2); hEff_long->SetLineWidth(2);  
97   hEff_long->SetXTitle("Pt, GeV/c"); 
98   hEff_long->SetYTitle("Efficiency"); 
99   
100   TH1F *hEff_short = new TH1F("hEff_short","Efficiency short",nPtSteps,minPt,maxPt); 
101   hEff_short->SetLineColor(8); hEff_short->SetLineWidth(2);  
102   hEff_short->SetXTitle("Pt, GeV/c"); 
103   hEff_short->SetYTitle("Efficiency"); 
104   
105   TH1F *hEff_long_eta = new TH1F("hEff_long_eta","Efficiency vs Eta",nEtaSteps,minEta,maxEta); 
106   hEff_long_eta->SetLineColor(2); hEff_long_eta->SetLineWidth(2);  
107   hEff_long_eta->SetXTitle("Pt, GeV/c"); 
108   hEff_long_eta->SetYTitle("Efficiency"); 
109   
110   TH1F *hEff_short_eta = new TH1F("hEff_short_eta","Efficiency short",nEtaSteps,minEta,maxEta); 
111   hEff_short_eta->SetLineColor(8); hEff_short_eta->SetLineWidth(2);  
112   hEff_short_eta->SetXTitle("Pt, GeV/c"); 
113   hEff_short_eta->SetYTitle("Efficiency"); 
114   
115   // page[4]
116   TH1F *hY=new TH1F("hY","Y resolution",50,-0.5,0.5);hY->SetLineColor(4);
117   hY->SetLineWidth(2);
118   hY->SetXTitle("(cm)");
119   TH1F *hZ=new TH1F("hZ","Z resolution",80,-4,4);hZ->SetLineColor(4);
120   hZ->SetLineWidth(2);
121   hZ->SetXTitle("(cm)");
122   TH1F *hPhi=new TH1F("hPhi","PHI resolution",100,-20.,20.); 
123   hPhi->SetFillColor(4);
124   hPhi->SetXTitle("(mrad)");
125   TH1F *hLambda=new TH1F("hLambda","Lambda resolution",100,-100,100);
126   hLambda->SetFillColor(17);
127   hLambda->SetXTitle("(mrad)");
128   TH1F *hPt=new TH1F("hPt","Relative Pt resolution",30,-10.,10.);  
129   
130   // page[5]
131   TH1F *hY_n=new TH1F("hY_n","Normalized Y resolution",50,-0.5,0.5); hY_n->SetLineColor(4);
132   hY_n->SetLineWidth(2);
133   hY_n->SetXTitle("  ");
134   TH1F *hZ_n=new TH1F("hZ_n","Normalized Z resolution",50,-0.5,0.5); hZ_n->SetLineColor(2);
135   hZ_n->SetLineWidth(2);
136   hZ_n->SetXTitle("   ");
137   TH1F *hLambda_n=new TH1F("hLambda_n","Normalized Lambda resolution",50,-0.5,0.5);
138   hLambda_n->SetFillColor(17);
139   TH1F *hPt_n=new TH1F("hPt_n","Normalized Pt resolution",50,-0.5,0.5); 
140   
141
142   Int_t nEvent = 0;
143
144   // Load Seeds saved as MC tracks 
145   TObjArray mctarray(2000);
146   TFile *mctf=TFile::Open("AliTRDtrackableSeeds.root");
147   if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDtrackableSeeds.root !\n"; return;}
148   TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
149   TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
150   Int_t const nMCtracks = (Int_t) mctracktree->GetEntries();
151   cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
152   for (Int_t i=0; i<nMCtracks; i++) {
153     AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack();
154     mctbranch->SetAddress(&ioMCtrack);
155     mctracktree->GetEvent(i);
156     mctarray.AddLast(ioMCtrack);
157   }
158   mctf->Close();                 
159
160   TFile *tf=TFile::Open("AliTRDtracks.root");
161
162   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
163   TObjArray rtarray(2000);
164
165   char   tname[100];
166   sprintf(tname,"TRDb_%d",nEvent);     
167   TTree *tracktree=(TTree*)tf->Get(tname);
168
169   TBranch *tbranch=tracktree->GetBranch("tracks");
170
171   Int_t nRecTracks = (Int_t) tracktree->GetEntries();
172   cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
173
174   for (Int_t i=0; i<nRecTracks; i++) {
175     AliTRDtrack *iotrack=new AliTRDtrack();
176     tbranch->SetAddress(&iotrack);
177     tracktree->GetEvent(i);
178     rtarray.AddLast(iotrack);
179   }
180   tf->Close();                 
181
182   //  return;
183
184   AliTRDmcTrack *mct;
185   AliTRDtrack *rt;
186   Int_t rtIndex[nMCtracks];
187   for(Int_t j = 0; j < nMCtracks; j++) {
188     mct = (AliTRDmcTrack*) mctarray.UncheckedAt(j);
189     rtIndex[j] = 0;
190     for (Int_t i=0; i<nRecTracks; i++) {
191       rt = (AliTRDtrack*) rtarray.UncheckedAt(i);
192       Int_t label = rt->GetSeedLabel();
193       if(mct->GetTrackIndex() == label) rtIndex[j] = i;
194     }
195   } 
196
197   // Load clusters
198
199   TFile *geofile =TFile::Open("AliTRDclusters.root");   
200   AliTRDtracker *Tracker = new AliTRDtracker(geofile);
201   Tracker->SetEventNumber(nEvent);
202
203   AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
204   AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
205   
206   Char_t *alifile = "AliTRDclusters.root";
207   TObjArray carray(2000);
208   TObjArray *ClustersArray = &carray;
209   Tracker->ReadClusters(ClustersArray,alifile);   
210
211
212   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
213   alifile = "galice.root";
214   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
215   if (!gafl) {
216     cout << "Open the ALIROOT-file " << alifile << endl;
217     gafl = new TFile(alifile);
218   }
219   else {
220     cout << alifile << " is already open" << endl;
221   }
222   gAlice = (AliRun*) gafl->Get("gAlice");
223   if (gAlice)
224     cout << "AliRun object found on file" << endl;
225   else
226     gAlice = new AliRun("gAlice","Alice test program");
227
228   AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");    
229   const Int_t nparticles = gAlice->GetEvent(nEvent);
230   if (nparticles <= 0) return;
231
232   TParticle *p;
233   Bool_t electrons[300000];
234   for(Int_t i = 0; i < 300000; i++) electrons[i] = kFALSE;
235
236   Bool_t mark_electrons = kFALSE;
237   if(mark_electrons) {
238     printf("mark electrons\n");
239
240     for(Int_t i = nPrimaries; i < nparticles; i++) {
241       p = gAlice->Particle(i);
242       if(p->GetMass() > 0.01) continue;
243       if(p->GetMass() < 0.00001) continue;
244       electrons[i] = kTRUE;
245     }
246   }
247
248   AliTRDcluster *cl = 0;
249   Int_t nw, label, index, ti[3];
250   Int_t mct_tbwc, rt_tbwc, mct_max_gap, rt_max_gap, mct_sector, rt_sector; 
251   Int_t mct_max_tb, rt_max_tb, mct_min_tb, rt_min_tb;
252
253   Int_t det, plane, ltb, gtb, gap, max_gap, sector;
254   Double_t Pt, Px, Py, Pz;
255   Double_t mcPt, mcPx, mcPy, mcPz;
256   Double_t x,y,z, mcX;
257   Int_t rtClusters, rtCorrect;
258
259   Int_t nToFind_long, nFound_long, nToFind_short, nFound_short;
260
261   printf("\n");
262
263   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
264   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
265   Double_t dx = (Double_t) fPar->GetTimeBinSize();   
266
267   Int_t tbAmp = fPar->GetTimeBefore();
268   Int_t maxAmp = 0;
269   Int_t tbDrift = fPar->GetTimeMax();
270   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
271
272   tbDrift = TMath::Min(tbDrift,maxDrift);
273   tbAmp = TMath::Min(tbAmp,maxAmp);
274
275   const Int_t nPlanes = fGeom->Nplan();
276   const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
277   const Int_t nTB = tbpp * nPlanes;
278   Int_t mask[nTB][2];
279
280   nToFind_long = 0; nFound_long = 0;
281   nToFind_short = 0; nFound_short = 0;
282
283   for(Int_t i=0; i < nMCtracks; i++) {
284     mct = (AliTRDmcTrack *) mctarray.UncheckedAt(i);
285     label = mct->GetTrackIndex();
286
287     // Waveform of the MC track
288     for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
289     Int_t mct_ncl = mct->GetNumberOfClusters();
290
291     for(ltb = 0; ltb < kMAX_TB; ltb++) {
292       for(plane = 0; plane < nPlanes; plane++) {
293         for(Int_t n = 0; n < 2; n++) {
294           index = mct->GetClusterIndex(ltb,plane,n);
295           if(index < 0) continue; 
296           cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
297
298           for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
299
300           if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
301             printf("mc wrong track label: %d, %d, %d != %d \n",
302                    ti[0], ti[1], ti[2], label);
303           } 
304           det=cl->GetDetector();
305           if(fGeom->GetPlane(det) != plane) 
306             printf("cluster plane = %d != %d expected plane\n", 
307                    fGeom->GetPlane(det), plane);
308           if(cl->GetLocalTimeBin() != ltb) 
309             printf("cluster ltb = %d != %d expected ltb\n", 
310                    cl->GetLocalTimeBin(), ltb);
311           gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
312           mask[gtb][n] = index;
313         }
314       }
315     }
316     
317     for(plane = 0; plane < nPlanes; plane++) {
318       for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
319         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
320         if(mask[gtb][0] > -1) break;
321       }
322       if(ltb >= -tbAmp) break;
323     }  
324     if((plane == nPlanes) && (ltb == -tbAmp-1)) {
325       // printf("warning: for track %d min tb is not found and set to %d!\n",
326       //         label, nTB-1);
327       mct_min_tb = nTB-1;
328       // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
329     }
330     else {
331       mct_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
332     }
333
334     for(plane = nPlanes-1 ; plane>=0; plane--) {
335       for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
336         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
337         if(mask[gtb][0] > -1) break;
338       }
339       if(ltb < tbDrift) break;
340     }  
341     if((plane == -1) && (ltb == tbDrift)) {
342       //      printf("warning: for track %d max tb is not found and set to 0!\n",label);
343       //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
344       mct_max_tb = 0;
345       //      mcX = Tracker->GetX(0,0,tbDrift-1);
346     }
347     else {
348       mct_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
349       //      mcX = Tracker->GetX(0,plane,ltb);
350       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
351       det=cl->GetDetector();
352       sector = fGeom->GetSector(det);  
353       mct_sector = sector;
354     }
355
356     mct_tbwc = 0;
357     mct_max_gap = 0;
358     gap = -1;
359     
360     for(nw = mct_min_tb; nw < mct_max_tb+1; nw++) {
361       gap++;
362       if(mask[nw][0] > -1) {
363         mct_tbwc++;
364         if(gap > mct_max_gap) mct_max_gap = gap;
365         gap = 0;
366       }
367     }
368
369     //  Waveform of the reconstructed track
370     if(rtIndex[i] >= 0) rt = (AliTRDtrack *) rtarray.UncheckedAt(rtIndex[i]);
371
372     for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
373     Int_t rt_ncl = rt->GetNumberOfClusters();
374
375     for(Int_t n = 0; n < rt_ncl; n++) {
376       index = rt->GetClusterIndex(n);
377       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
378       
379       for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
380       
381       if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
382         // printf("rt wrong track label: %d, %d, %d != %d \n", ti[0], ti[1], ti[2], label);
383         continue;
384       } 
385       
386       det=cl->GetDetector();
387       plane = fGeom->GetPlane(det); 
388       ltb = cl->GetLocalTimeBin(); 
389       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
390       mask[gtb][0] = index;
391     }
392
393     for(plane = 0; plane < nPlanes; plane++) {
394       for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
395         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
396         if(mask[gtb][0] > -1) break;
397       }
398       if(ltb >= -tbAmp) break;
399     }  
400     if((plane == nPlanes) && (ltb == -tbAmp-1)) {
401       // printf("warning: for track %d min tb is not found and set to %d!\n",
402       //             label, nTB-1);
403       rt_min_tb = nTB-1;
404       // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
405     }
406     else {
407       rt_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
408     }
409
410     for(plane = nPlanes-1 ; plane>=0; plane--) {
411       for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
412         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
413         if(mask[gtb][0] > -1) break;
414       }
415       if(ltb < tbDrift) break;
416     }  
417     if((plane == -1) && (ltb == tbDrift)) {
418       // printf("warning: for track %d max tb is not found and set to 0!\n",label);
419       //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
420       rt_max_tb = 0;
421       //      mcX = Tracker->GetX(0,0,tbDrift-1);
422     }
423     else {
424       rt_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
425       //      mcX = Tracker->GetX(0,plane,ltb);
426       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
427       det=cl->GetDetector();
428       sector = fGeom->GetSector(det);  
429       rt_sector = sector;
430     }
431
432     rt_tbwc = 0;
433     rt_max_gap = 0;
434     gap = -1;
435     
436     for(nw = rt_min_tb; nw < rt_max_tb+1; nw++) {
437       gap++;
438       if(mask[nw][0] > -1) {
439         rt_tbwc++;
440         if(gap > rt_max_gap) rt_max_gap = gap;
441         gap = 0;
442       }
443     }
444
445     // Fill the histoes
446
447     if(page[0]) {
448       hNcl->Fill((Float_t) mct_ncl);
449       hNtb->Fill((Float_t) mct_tbwc);
450       hNep->Fill((Float_t) mct_max_tb);
451       hNmg->Fill((Float_t) mct_max_gap);
452     }
453     if(page[1]) {
454       h2ep->Fill((Float_t) rt_max_tb, (Float_t) mct_max_tb);
455       h2ntb->Fill((Float_t) rt_tbwc, (Float_t) mct_tbwc);
456       h2mg->Fill((Float_t) rt_max_gap, (Float_t) mct_max_gap);
457     }
458     if(page[2]) {
459       p = gAlice->Particle(label);
460       hPt_all->Fill(p->Pt());
461       hEta_all->Fill(p->Eta());
462       if(mct_max_tb > 60) {
463         nToFind_long++;
464         hPt_long->Fill(p->Pt());
465         hEta_long->Fill(p->Eta());
466         if(((mct_max_tb - rt_max_tb) < 10) && 
467            (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
468           nFound_long++;
469           hrtPt_long->Fill(p->Pt());
470           hrtEta_long->Fill(p->Eta());
471         }         
472       }
473       if((mct_max_tb < 60) && (mct_max_tb > 10)) {
474         nToFind_short++;
475         hPt_short->Fill(p->Pt());
476         hEta_short->Fill(p->Eta());
477         if(((mct_max_tb - rt_max_tb) < 10) && 
478            (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
479           nFound_short++;
480           hrtPt_short->Fill(p->Pt());
481           hrtEta_short->Fill(p->Eta());
482         }         
483       }
484     }
485     if(page[4] && page[5]) {
486       if((mct_tbwc > 50) && (rt_tbwc > 50)) {
487         if(rt->GetSeedLabel() != label) {
488           printf("mct vs rt index mismatch: %d != %d \n",
489                  label, rt->GetSeedLabel());
490           return;
491         }
492         Pt = rt->GetPt(); 
493         Double_t cc = TMath::Abs(rt->GetSigmaC2()); 
494         mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
495         rt->GetPxPyPz(Px,Py,Pz);      
496
497         printf("\n\ntrack %d \n", label);
498         printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); 
499         printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); 
500     
501         mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
502         if(mcPt < 0.0001) mcPt = 0.0001;
503     
504         Float_t lamg=TMath::ATan2(mcPz,mcPt);
505         Float_t lam =TMath::ATan2(Pz,Pt);
506         if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
507         else hPt_n->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
508     
509         Float_t phig=TMath::ATan2(mcPy,mcPx);
510         Float_t phi =TMath::ATan2(Py,Px);    
511         //      if(!(rt->PropagateTo(x))) continue;
512         //    if((mcPt > Pt_min) && (mcPt < Pt_max)) {
513         hLambda->Fill((lam - lamg)*1000.);
514         hLambda_n->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
515         if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
516         else hPt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); 
517         hPhi->Fill((phi - phig)*1000.);
518         hY->Fill(rt->GetY() - y);
519         hZ->Fill(rt->GetZ() - z);
520         hY_n->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
521         hZ_n->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
522       }
523     }  
524   }
525
526   // plot pictures
527
528   if(page[0]) {
529     TCanvas *c0=new TCanvas("c0","",0,0,700,850);
530     gStyle->SetOptStat(111110);
531     c0->Divide(2,2);
532     c0->cd(1); gPad->SetLogy(); hNcl->Draw();
533     c0->cd(2); hNtb->Draw();
534     c0->cd(3); hNep->Draw();
535     c0->cd(4); hNmg->Draw();
536     c0->Print("c0.ps","ps"); 
537   }
538   if(page[1]) {
539     TCanvas *c1=new TCanvas("c1","",0,0,700,850);
540     gStyle->SetOptStat(0);
541     c1->Divide(2,2);
542     c1->cd(1); h2ep->Draw();
543     c1->cd(2); h2ntb->Draw();
544     c1->cd(3); h2mg->Draw();
545     //    c1->cd(4); hNmg->Draw();
546     c1->Print("c1.ps","ps"); 
547   }
548   if(page[2]) {
549     TCanvas *c2=new TCanvas("c2","",0,0,700,850);
550     gStyle->SetOptStat(0);
551     c2->Divide(2,2);
552     c2->cd(1); hPt_all->Draw(); hPt_long->Draw("same"); hPt_short->Draw("same"); 
553     c2->cd(2); hEta_all->Draw(); hEta_long->Draw("same"); hEta_short->Draw("same"); 
554     //    c2->cd(3); h2mg->Draw();
555     //    c2->cd(4); hNmg->Draw();
556     c2->Print("c2.ps","ps"); 
557   }
558   if(page[3]) {
559     TCanvas *c3=new TCanvas("c3","",0,0,700,850);
560     gStyle->SetOptStat(0);
561     c3->Divide(1,2);
562     c3->cd(1);
563     hrtPt_long->Sumw2(); hPt_long->Sumw2(); 
564     hEff_long->Divide(hrtPt_long,hPt_long,1,1.,"b");
565     hEff_long->SetMaximum(1.4);
566     hEff_long->SetYTitle("Matching Efficiency");
567     hEff_long->SetXTitle("Pt (GeV/c)");
568     hEff_long->Draw();          
569     hrtPt_short->Sumw2(); hPt_short->Sumw2(); 
570     hEff_short->Divide(hrtPt_short,hPt_short,1,1.,"b");
571     hEff_short->SetMaximum(1.4);
572     hEff_short->SetYTitle("Matching Efficiency");
573     hEff_short->SetXTitle("Pt (GeV/c)");
574     hEff_short->Draw("same");          
575
576     c3->cd(2);
577     hrtEta_long->Sumw2(); hEta_long->Sumw2(); 
578     hEff_long_eta->Divide(hrtEta_long,hEta_long,1,1.,"b");
579     hEff_long_eta->SetMaximum(1.4);
580     hEff_long_eta->SetYTitle("Matching Efficiency");
581     hEff_long_eta->SetXTitle("Eta");
582     hEff_long_eta->Draw();          
583     hrtEta_short->Sumw2(); hEta_short->Sumw2(); 
584     hEff_short_eta->Divide(hrtEta_short,hEta_short,1,1.,"b");
585     hEff_short_eta->SetMaximum(1.4);
586     hEff_short_eta->SetYTitle("Matching Efficiency");
587     hEff_short_eta->SetXTitle("Eta");
588     hEff_short_eta->Draw("same");          
589     c3->Print("c3.ps","ps"); 
590   }
591   if(page[4]) {
592     TCanvas *c4=new TCanvas("c4","",0,0,700,850);
593     gStyle->SetOptStat(111110);
594     c4->Divide(2,3);
595     c4->cd(1); hY->Draw();
596     c4->cd(2); hZ->Draw();
597     c4->cd(3); hPhi->Draw();
598     c4->cd(4); hLambda->Draw();
599     c4->cd(5); hPt->Draw();
600     c4->Print("c4.ps","ps"); 
601   }
602   if(page[5]) {
603     TCanvas *c5=new TCanvas("c5","",0,0,700,850);
604     gStyle->SetOptStat(111110);
605     c5->Divide(2,3);
606     c5->cd(1); hY_n->Draw();
607     c5->cd(2); hZ_n->Draw();
608     c5->cd(3); hLambda_n->Draw();
609     c5->cd(4); hPt_n->Draw();
610     c5->Print("c5.ps","ps"); 
611   }
612   printf("Efficiency (long) = %d/%d = %f\n",nFound_long,nToFind_long,
613          ((Float_t) nFound_long / (Float_t) nToFind_long));
614   printf("Efficiency (shrt) = %d/%d = %f\n",nFound_short,nToFind_short,
615          ((Float_t) nFound_short / (Float_t) nToFind_short));
616 }
617
618