TRD/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 = TMath::Abs(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