TRD/AliTRDbackTrackAnalysis.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 void AliTRDbackTrackAnalysis() {
  18
  19   //  const Int_t nPrimaries = 84210/400;
  20   const Int_t nPrimaries = 84210/16;
  21
  22
  23   Float_t Pt_min = 0.;
  24   Float_t Pt_max = 20.;
  25
  26   TH1F *hp=new TH1F("hp","PHI resolution",100,-20.,20.); hp->SetFillColor(4);
  27   TH1F *hl=new TH1F("hl","LAMBDA resolution",100,-100,100);hl->SetFillColor(4);
  28   TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
  29
  30   TH1F *hpta=new TH1F("hpta","norm. Pt resolution",50,-5.,5.);
  31   hpta->SetFillColor(17);
  32   hpta->SetXTitle("(%)");
  33
  34   TH1F *hla=new TH1F("hla","norm. Lambda resolution",50,-5.,5.);
  35   hla->SetFillColor(17);
  36   TH1F *hya=new TH1F("hya","norm. Y resolution",50,-5.,5.);
  37   hya->SetFillColor(17);
  38   TH1F *hza=new TH1F("hza","norm. Z resolution",50,-5.,5.);
  39   hza->SetFillColor(17);
  40
  41   hpt->SetFillColor(2);
  42
  43   TH1F *hy=new TH1F("hy","Y resolution",50,-0.5,0.5);hy->SetLineColor(4);
  44   hy->SetLineWidth(2);
  45   hy->SetXTitle("(cm)");
  46
  47   TH1F *hz=new TH1F("hz","Z resolution",80,-4,4);hz->SetLineColor(2);
  48   hz->SetLineWidth(2);
  49   hz->SetXTitle("(cm)");
  50
  51   TH1F *hx=new TH1F("hx","X(out)",150,250.,400.); hx->SetFillColor(4);
  52
  53   const Int_t nPtSteps = 30;
  54   const Float_t maxPt = 3.;
  55
  56   TH1F *hgood=new TH1F("hgood","long TRD tracks, TPC seeds",nPtSteps,0,maxPt);
  57   hgood->SetYTitle("Counts");
  58   hgood->SetXTitle("Pt (GeV/c)");
  59   hgood->SetLineColor(3);
  60
  61   TH1F *hGoodAndSeed = new TH1F("GoodAndSeed","TPC seed and good",nPtSteps,0,maxPt);
  62   hGoodAndSeed->SetLineColor(2);
  63
  64   TH2F *hRTvsMC = new TH2F("RTvsMC","RTvsMC",100,-4.5,95.5,100,-4.5,95.5);
  65   hRTvsMC->SetMarkerColor(4);
  66   hRTvsMC->SetMarkerSize(2);
  67
  68   TH2F *hXvsMCX = new TH2F("XvsMCX","XvsMCX",150,250,400,150,250,400);
  69   hXvsMCX->SetMarkerColor(4);
  70   hXvsMCX->SetMarkerSize(2);
  71
  72   TH2F *h2seed = new TH2F("2dGood","TPC seeds",60,0.,60,50,0.,3.);
  73   h2seed->SetMarkerColor(4);
  74   h2seed->SetMarkerSize(2);
  75
  76   TH2F *h2lost = new TH2F("2dSeedAndGood","SeedButNotGood",60,0.,60.,50,0.,3.);
  77   h2lost->SetMarkerColor(2);
  78   h2lost->SetMarkerSize(2);
  79
  80
  81   TH1F *hseed=new TH1F("seed","TPC seeds",nPtSteps,0,maxPt);
  82   hseed->SetLineColor(4);
  83
  84
  85   TH1F *hfound=new TH1F("hfound","Found tracks",nPtSteps,0,maxPt);
  86   hfound->SetYTitle("Counts");
  87   hfound->SetXTitle("Pt (GeV/c)");
  88
  89   TH1F *heff=new TH1F("heff","Matching Efficiency",nPtSteps,0,maxPt); // efficiency, good tracks
  90   heff->SetLineColor(4); heff->SetLineWidth(2);
  91   heff->SetXTitle("Pt, GeV/c");
  92   heff->SetYTitle("Efficiency");
  93
  94   TH1F *hSeedEff = new TH1F("hSeedEff","TPC Efficiency",nPtSteps,0,maxPt);
  95   hSeedEff->SetLineColor(4); hSeedEff->SetLineWidth(2);
  96   hSeedEff->SetXTitle("Pt, GeV/c");
  97   hSeedEff->SetYTitle("Efficiency");
  98
  99
 100   TH1F *hol=new TH1F("hol","Overlap fraction",105,-2.5,102.5);
 101   hol->SetLineColor(4); hol->SetLineWidth(2);
 102   hol->SetXTitle("Fraction,(%)");
 103   hol->SetYTitle("Counts");
 104
 105   TH1F *hend=new TH1F("end","missing tail",80,-10.5,69.5);
 106
 107   TH1F *hFraction=new TH1F("fraction","Fraction of found clusters",110,0,1.1);
 108   TH1F *hCorrect=new TH1F("correct","Fraction of correct clusters",110,0,1.1);
 109
 110
 111   Int_t nEvent = 0;
 112   const Int_t maxIndex = nPrimaries;
 113   Bool_t seedLabel[maxIndex];
 114   Int_t mcIndex[maxIndex];
 115   Int_t rtIndex[maxIndex];
 116
 117   for(Int_t i = 0; i < maxIndex; i++) {
 118     seedLabel[i] = kFALSE;
 119     mcIndex[i] = -1;
 120     rtIndex[i] = -1;
 121   }
 122
 123   // mark available seeds from TPC
 124   Int_t nSeeds = 0;
 125   Int_t nPrimarySeeds = 0;
 126
 127   printf("marking found seeds from TPC\n");
 128   TDirectory *savedir=gDirectory;
 129
 130   TFile *in=TFile::Open("AliTPCBackTracks.root");
 131   if (!in->IsOpen()) {
 132     cerr<<"can't open file AliTPCBackTracks.root  !\n"; return;
 133   }
 134
 135   char   tname[100];
 136   sprintf(tname,"seedsTPCtoTRD_%d",nEvent);
 137   TTree *seedTree=(TTree*)in->Get(tname);
 138   if (!seedTree) {
 139      cerr<<"AliTRDtracker::PropagateBack(): ";
 140      cerr<<"can't get a tree with seeds from TPC !\n";
 141   }
 142
 143   AliTPCtrack *seed=new AliTPCtrack;
 144   seedTree->SetBranchAddress("tracks",&seed);
 145
 146   Int_t n=(Int_t)seedTree->GetEntries();
 147   for (Int_t i=0; i<n; i++) {
 148      seedTree->GetEvent(i);
 149      Int_t lbl = seed->GetLabel();
 150      if(lbl < 0) {
 151        printf("negative seed label %d \n",lbl);
 152        continue;
 153      }
 154      if(lbl >= maxIndex) continue;
 155      seedLabel[lbl] = kTRUE;
 156      if(lbl < nPrimaries) nPrimarySeeds++;
 157      nSeeds++;
 158   }
 159   delete seed;
 160   delete seedTree;
 161
 162   printf("Found %d seeds from primaries among overall %d seeds \n",
 163          nPrimarySeeds, nSeeds);
 164
 165   savedir->cd();
 166   // done with marking TPC seeds
 167
 168
 169   TFile *tf=TFile::Open("AliTRDtracks.root");
 170
 171   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
 172   TObjArray tarray(2000);
 173
 174   sprintf(tname,"TRDb_%d",nEvent);
 175   TTree *tracktree=(TTree*)tf->Get(tname);
 176
 177   TBranch *tbranch=tracktree->GetBranch("tracks");
 178
 179   Int_t nRecTracks = (Int_t) tracktree->GetEntries();
 180   cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
 181
 182   for (Int_t i=0; i<nRecTracks; i++) {
 183     AliTRDtrack *iotrack=new AliTRDtrack();
 184     tbranch->SetAddress(&iotrack);
 185     tracktree->GetEvent(i);
 186     tarray.AddLast(iotrack);
 187     Int_t trackLabel = iotrack->GetLabel();
 188
 189     //    printf("rt with %d clusters and label %d \n",
 190     //     iotrack->GetNumberOfClusters(), trackLabel);
 191
 192     if(trackLabel < 0) continue;
 193     if(trackLabel >= maxIndex) continue;
 194     rtIndex[trackLabel] = i;
 195   }
 196   tf->Close();
 197
 198   //  return;
 199
 200   // Load MC tracks
 201   TFile *mctf=TFile::Open("AliTRDmcTracks.root");
 202   if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDmcTracks.root !\n"; return;}
 203   TObjArray mctarray(2000);
 204   TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
 205   TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
 206   Int_t nMCtracks = (Int_t) mctracktree->GetEntries();
 207   cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
 208   for (Int_t i=0; i<nMCtracks; i++) {
 209     AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack;
 210     mctbranch->SetAddress(&ioMCtrack);
 211     mctracktree->GetEvent(i);
 212     mctarray.AddLast(ioMCtrack);
 213     Int_t mcLabel = ioMCtrack->GetTrackIndex();
 214     if(mcLabel < 0) {printf("negative mc label detected!\n"); continue;}
 215     if(mcLabel >= maxIndex) continue;
 216     mcIndex[mcLabel] = i;
 217   }
 218   mctf->Close();
 219
 220
 221   // Load clusters
 222
 223   TFile *geofile =TFile::Open("AliTRDclusters.root");
 224   AliTRDtracker *Tracker = new AliTRDtracker(geofile);
 225   Tracker->SetEventNumber(nEvent);
 226
 227   AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry");
 228   AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter");
 229
 230   Char_t *alifile = "AliTRDclusters.root";
 231   TObjArray carray(2000);
 232   TObjArray *ClustersArray = &carray;
 233   Tracker->ReadClusters(ClustersArray,alifile);
 234
 235
 236   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
 237   alifile = "galice.root";
 238
 239   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
 240   if (!gafl) {
 241     cout << "Open the ALIROOT-file " << alifile << endl;
 242     gafl = new TFile(alifile);
 243   }
 244   else {
 245     cout << alifile << " is already open" << endl;
 246   }
 247
 248   // Get AliRun object from file or create it if not on file
 249   gAlice = (AliRun*) gafl->Get("gAlice");
 250   if (gAlice)
 251     cout << "AliRun object found on file" << endl;
 252   else
 253     gAlice = new AliRun("gAlice","Alice test program");
 254
 255
 256   // Define the objects
 257   AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");
 258
 259   // Import the Trees for the event nEvent in the file
 260   const Int_t nparticles = gAlice->GetEvent(nEvent);
 261   if (nparticles <= 0) return;
 262
 263   TParticle *p;
 264   Bool_t electrons[300000] = { kFALSE };
 265
 266   Bool_t mark_electrons = kFALSE;
 267   if(mark_electrons) {
 268     printf("mark electrons\n");
 269
 270     for(Int_t i = nPrimaries; i < nparticles; i++) {
 271       p = gAlice->Particle(i);
 272       if(p->GetMass() > 0.01) continue;
 273       if(p->GetMass() < 0.00001) continue;
 274       electrons[i] = kTRUE;
 275     }
 276   }
 277
 278   AliTRDcluster *cl = 0;
 279   Int_t nw, label, index, ti[3], tbwc;
 280   Int_t det, plane, ltb, gtb, gap, max_gap, sector, mc_sector, min_tb, max_tb;
 281   Double_t Pt, Px, Py, Pz;
 282   Double_t mcPt, mcPx, mcPy, mcPz;
 283   Double_t x,y,z, mcX;
 284   Int_t rtClusters, rtCorrect;
 285
 286   printf("\n");
 287
 288   AliTRDmcTrack *mct = 0;
 289   AliTRDtrack *rt = 0;
 290
 291   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
 292   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
 293   Double_t dx = (Double_t) fPar->GetTimeBinSize();
 294
 295   Int_t tbAmp = fPar->GetTimeBefore();
 296   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
 297   Int_t tbDrift = fPar->GetTimeMax();
 298   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
 299
 300   tbDrift = TMath::Min(tbDrift,maxDrift);
 301   tbAmp = TMath::Min(tbAmp,maxAmp);
 302
 303   const Int_t nPlanes = fGeom->Nplan();
 304   const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
 305   const Int_t nTB = tbpp * nPlanes;
 306   Int_t mask[nTB];
 307
 308   for(Int_t i=0; i < maxIndex; i++) {
 309
 310     rt = 0; mct = 0;
 311     if(rtIndex[i] >= 0) rt = (AliTRDtrack *) tarray.UncheckedAt(rtIndex[i]);
 312     if(mcIndex[i] >= 0) mct = (AliTRDmcTrack *) mctarray.UncheckedAt(mcIndex[i]);
 313
 314     if(!mct) continue;
 315     label = mct->GetTrackIndex();
 316
 317     //    if(TMath::Abs(mct->GetMass()-0.136) > 0.01) continue;
 318
 319     Int_t ncl = mct->GetNumberOfClusters();
 320
 321     // check how many time bins have a cluster
 322
 323     for(nw = 0; nw < nTB; nw++) mask[nw] = -1;
 324
 325     for(Int_t j = 0; j < ncl; j++) {
 326       index = mct->GetClusterIndex(j);
 327       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
 328
 329       for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
 330
 331       if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
 332         printf("wrong track label: %d, %d, %d != %d \n",
 333                ti[0], ti[1], ti[2], label);
 334       }
 335       det=cl->GetDetector();
 336       plane = fGeom->GetPlane(det);
 337       ltb = cl->GetLocalTimeBin();
 338       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
 339       mask[gtb] = index;
 340     }
 341
 342
 343
 344
 345     for(plane = 0; plane < nPlanes; plane++) {
 346       for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
 347         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
 348         if(mask[gtb] > -1) break;
 349       }
 350       if(ltb >= -tbAmp) break;
 351     }
 352     if((plane == nPlanes) && (ltb == -tbAmp-1)) {
 353       printf("warning: for track %d min tb is not found and set to %d!\n",
 354              label, nTB-1);
 355       min_tb = nTB-1;
 356       //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
 357     }
 358     else {
 359       min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
 360     }
 361
 362     for(plane = nPlanes-1 ; plane>=0; plane--) {
 363       for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
 364         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
 365         if(mask[gtb] > -1) break;
 366       }
 367       if(ltb < tbDrift) break;
 368     }
 369     if((plane == -1) && (ltb == tbDrift)) {
 370       printf("warning: for track %d max tb is not found and set to 0!\n",label);
 371       //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
 372       max_tb = 0;
 373       mcX = Tracker->GetX(0,0,tbDrift-1);
 374     }
 375     else {
 376       max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
 377       mcX = Tracker->GetX(0,plane,ltb);
 378       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb]);
 379       det=cl->GetDetector();
 380       sector = fGeom->GetSector(det);
 381       mc_sector = sector;
 382     }
 383
 384
 385     tbwc = 0;
 386     max_gap = 0;
 387     gap = -1;
 388
 389     for(nw = min_tb; nw < max_tb+1; nw++) {
 390       gap++;
 391       if(mask[nw] > -1) {
 392         tbwc++;
 393         if(gap > max_gap) max_gap = gap;
 394         gap = 0;
 395       }
 396     }
 397
 398     if(tbwc < ((Int_t) (nTB * Tracker->GetMinClustersInTrack()))) continue;
 399     if(gap > Tracker->GetMaxGap()) continue;
 400     p = gAlice->Particle(label);
 401
 402     printf("good track %d with min_tb %d, max_tb %d, gap %d\n",
 403            label,min_tb,max_tb,gap);
 404
 405     hgood->Fill(p->Pt());
 406
 407     if(!rt) continue;
 408
 409     if(rt->GetLabel() != label) {
 410       printf("mct vs rt index mismatch: %d != %d \n",
 411              label, rt->GetLabel());
 412       return;
 413     }
 414
 415     if(!seedLabel[i]) continue;
 416
 417     hGoodAndSeed->Fill(p->Pt());
 418
 419     hXvsMCX->Fill(mcX,rt->GetX());
 420
 421     rtClusters = rt->GetNumberOfClusters();
 422
 423     // find number of tb with correct cluster
 424     rtCorrect = 0;
 425
 426     for(Int_t j = 0; j < rtClusters; j++) {
 427       index = rt->GetClusterIndex(j);
 428       cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
 429
 430       for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
 431
 432       if((ti[0] != label)&&(ti[1] != label)&&(ti[2] != label)) continue;
 433       rtCorrect++;
 434     }
 435
 436     Float_t  foundFraction = ((Float_t) rtCorrect) / (tbwc + 0.00001);
 437     Float_t  correctFraction = ((Float_t) rtCorrect) / ((Float_t) rtClusters);
 438
 439     if(foundFraction > 1) printf("fraction = %f/%f \n",
 440                                  (Float_t) rtCorrect, (Float_t) tbwc);
 441
 442
 443     if(foundFraction <= 1) hFraction->Fill(foundFraction);
 444     hCorrect->Fill(correctFraction);
 445
 446     hRTvsMC->Fill((Float_t) tbwc, (Float_t) rtClusters);
 447
 448     if((foundFraction < 0.7) || (correctFraction < 0.7)) {
 449       printf("not found track %d with FrctnFound %f and FrctnCorrect %f\n",
 450              label, foundFraction, correctFraction);
 451       continue;
 452     }
 453
 454     hfound->Fill(p->Pt());
 455
 456     Pt = TMath::Abs(rt->GetPt());
 457     Double_t cc = TMath::Abs(rt->GetSigmaC2());
 458     mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
 459     rt->GetPxPyPz(Px,Py,Pz);
 460
 461     printf("\n\ntrack %d \n", label);
 462     printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz);
 463     printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz);
 464
 465     mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
 466     if(mcPt < 0.0001) mcPt = 0.0001;
 467
 468     Float_t lamg=TMath::ATan2(mcPz,mcPt);
 469     Float_t lam =TMath::ATan2(Pz,Pt);
 470     if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
 471     else hpta->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
 472
 473     Float_t phig=TMath::ATan2(mcPy,mcPx);
 474     Float_t phi =TMath::ATan2(Py,Px);
 475
 476
 477     if(!(rt->PropagateTo(x))) continue;
 478
 479
 480     if((mcPt > Pt_min) && (mcPt < Pt_max)) {
 481       hl->Fill((lam - lamg)*1000.);
 482       hla->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
 483       if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
 484       else hpt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.);
 485       hp->Fill((phi - phig)*1000.);
 486       hy->Fill(rt->GetY() - y);
 487       hz->Fill(rt->GetZ() - z);
 488       hya->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
 489       hza->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
 490       hx->Fill((Float_t) x);
 491     }
 492   }
 493
 494
 495   gStyle->SetOptStat(0);
 496   //  gStyle->SetOptStat(111110);
 497   gStyle->SetOptFit(1);
 498
 499
 500   /*
 501   TCanvas *c1=new TCanvas("c1","",0,0,700,850);
 502
 503   //  gStyle->SetOptStat(0);
 504   //  gStyle->SetOptStat(111110);
 505   gStyle->SetOptFit(1);
 506
 507   TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
 508   p1->cd(); p1->SetFillColor(10); p1->SetFrameFillColor(10);
 509   hgood->Draw(); hGoodAndSeed->Draw("same"); c1->cd();
 510
 511   TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
 512   p2->cd(); p2->SetFillColor(10); p2->SetFrameFillColor(10);
 513   hFraction->SetYTitle("Counts");
 514   hFraction->SetXTitle("Fraction");
 515   hFraction->Draw(); c1->cd();
 516
 517   TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
 518   p3->cd(); p3->SetFillColor(10); p3->SetFrameFillColor(10);
 519   hfound->Draw(); c1->cd();
 520
 521   TPad *p4=new TPad("p4","",0.5,0,1.,0.3); p4->Draw();
 522   p4->cd(); p4->SetFillColor(10); p4->SetFrameFillColor(10);
 523   hCorrect->SetYTitle("Counts");
 524   hCorrect->SetXTitle("Fraction");
 525   hCorrect->Draw(); c1->cd();
 526
 527   TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
 528   p5->SetFillColor(10); p5->SetFrameFillColor(10);
 529   hgood->Sumw2(); hGoodAndSeed->Sumw2(); // hfake->Sumw2();
 530   hSeedEff->Divide(hGoodAndSeed,hgood,1,1.,"b");
 531   //  hf->Divide(hfake,hgood,1,1.,"b");
 532   hSeedEff->SetMaximum(1.4);
 533   hSeedEff->SetYTitle("TPC efficiency");
 534   hSeedEff->SetXTitle("Pt (GeV/c)");
 535   hSeedEff->Draw();
 536
 537   TLine *line1 = new TLine(0,1.0,maxPt,1.0); line1->SetLineStyle(4);
 538   line1->Draw("same");
 539   TLine *line2 = new TLine(0,0.9,maxPt,0.9); line2->SetLineStyle(4);
 540   line2->Draw("same");
 541   */
 542
 543   TCanvas *c2=new TCanvas("c2","",0,0,700,850);
 544
 545   TPad *p12=new TPad("p12","",0,0.3,.5,.6); p12->Draw();
 546   p12->cd(); p12->SetFillColor(42); p12->SetFrameFillColor(10);
 547   hp->SetFillColor(4);  hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c2->cd();
 548
 549   TPad *p22=new TPad("p22","",0.5,.3,1,.6); p22->Draw();
 550   p22->cd(); p22->SetFillColor(42); p22->SetFrameFillColor(10);
 551   hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c2->cd();
 552
 553   TPad *p32=new TPad("p32","",0,0,0.5,0.3); p32->Draw();
 554   p32->cd(); p32->SetFillColor(42); p32->SetFrameFillColor(10);
 555   hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c2->cd();
 556
 557   TPad *p42=new TPad("p42","",0.5,0,1.,0.3); p42->Draw();
 558   p42->cd(); p42->SetFillColor(42); p42->SetFrameFillColor(10);
 559   hgood->Draw(); hGoodAndSeed->Draw("same"); c2->cd();
 560
 561   TPad *p52=new TPad("p52","",0,0.6,1,1); p52->Draw(); p52->cd();
 562   p52->SetFillColor(41); p52->SetFrameFillColor(10);
 563   hfound->Sumw2();
 564   heff->Divide(hfound,hGoodAndSeed,1,1.,"b");
 565   //  hf->Divide(hfake,hgood,1,1.,"b");
 566   heff->SetMaximum(1.4);
 567   heff->SetXTitle("Pt (GeV/c)");
 568   heff->Draw();
 569
 570   TLine *line12 = new TLine(0,1.0,maxPt,1.0); line12->SetLineStyle(4);
 571   line12->Draw("same");
 572   TLine *line22 = new TLine(0,0.9,maxPt,0.9); line22->SetLineStyle(4);
 573   line22->Draw("same");
 574
 575   c2->Print("matching.ps","ps");
 576
 577
 578   /*
 579   TCanvas *cxyz = new TCanvas("cxyz","",50,50,750,900);
 580   cxyz->Divide(2,2);
 581   cxyz->cd(1); hx->Draw();
 582   cxyz->cd(2); hy->Draw();
 583   cxyz->cd(3); hz->Draw();
 584   cxyz->cd(4); hXvsMCX->Draw();
 585   */
 586
 587   TCanvas *cs = new TCanvas("cs","",0,0,700,850);
 588   cs->Divide(2,3);
 589
 590   cs->cd(1); hy->Fit("gaus");
 591   cs->cd(2); hz->Fit("gaus");
 592   cs->cd(3); hpta->Fit("gaus");
 593   cs->cd(4); hla->Fit("gaus");
 594   cs->cd(5); hya->Fit("gaus");
 595   cs->cd(6); hza->Fit("gaus");
 596
 597   cs->Print("resolution.ps","ps");
 598
 599   /*
 600   TCanvas *cvs = new TCanvas("cvs","",0,0,700,850);
 601   cvs->cd(); hRTvsMC->Draw();
 602   */
 603
 604 }
 605
 606