PWG1/comparison/AliTPCComparisonTask.cxx

   1 //-------------------------------------------------------------------------
   2 //
   3 // This is the PROOF-enabled version of TPC/AliTPCComparison.C macro.
   4 // Origin:  Andrei.Zalite@cern.ch
   5 //
   6 //-------------------------------------------------------------------------
   7
   8 #include "TChain.h"
   9 #include "TTree.h"
  10 #include "TH1F.h"
  11 #include "TH2F.h"
  12 #include "TList.h"
  13 #include "TClonesArray.h"
  14 #include "TCanvas.h"
  15 #include "TStyle.h"
  16 #include "TLine.h"
  17 #include "TText.h"
  18 #include "TFile.h"
  19
  20 #include "AliLog.h"
  21 #include "AliVEvent.h"
  22 #include "AliESDEvent.h"
  23 #include "AliMCEvent.h"
  24 #include "AliESDtrack.h"
  25 #include "AliTrackReference.h"
  26 #include "AliMCComparisonTrack.h"
  27 #include "AliMCComparisonTrack.cxx"
  28
  29 #include "AliTPCComparisonTask.h"
  30
  31 ClassImp(AliTPCComparisonTask)
  32
  33
  34 AliTPCComparisonTask::AliTPCComparisonTask()
  35   : AliAnalysisTaskSE("AliTPCComparisonTask"),
  36     fListOfHistos(0),
  37     fGood(0),
  38     fFound(0),
  39     fFake(0),
  40     fP(0),
  41     fL(0),
  42     fPt(0),
  43     fHmpt(0),
  44     fE(0),
  45     fEp(0),
  46     fGoodPhi(0),
  47     fFoundPhi(0)
  48 {
  49   // Default constructor
  50   AliInfo("Default constructor AliTPCComparisonTask");
  51   // Define input and output slots here
  52   // Input slot #0 works with a TChain
  53   DefineInput(0, TChain::Class());
  54   // Output slot #1 TList
  55   DefineOutput(1, TList::Class());
  56 }
  57
  58
  59 AliTPCComparisonTask::AliTPCComparisonTask(const char* name)
  60   : AliAnalysisTaskSE(name),
  61     fListOfHistos(0),
  62     fGood(0),
  63     fFound(0),
  64     fFake(0),
  65     fP(0),
  66     fL(0),
  67     fPt(0),
  68     fHmpt(0),
  69     fE(0),
  70     fEp(0),
  71     fGoodPhi(0),
  72     fFoundPhi(0)
  73 {
  74   // Constructor
  75   AliInfo("Constructor AliTPCComparisonTask");
  76   // Define input and output slots here
  77   // Input slot #0 works with a TChain
  78   DefineInput(0, TChain::Class());
  79   // Output slot #1 TList
  80   DefineOutput(1, TList::Class());
  81 }
  82
  83
  84 void AliTPCComparisonTask::UserCreateOutputObjects()
  85 {
  86   // Create histograms
  87   // Called once
  88   AliInfo("AliTPCComparisonTask::UserCreateOutputObjects");
  89   // Create output container
  90   fListOfHistos = new TList();
  91
  92   fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
  93   fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
  94   fFake = new TH1F("fFake",  "Pt for fake tracks", 34, 0.2, 7.0);
  95   fP = new TH1F("fP", "PHI resolution", 50, -20., 20.);
  96   fL = new TH1F("fL", "LAMBDA resolution", 50, -20., 20.);
  97   fPt = new TH1F("fPt", "Relative Pt resolution", 30, -10., 10.);
  98   fHmpt = new TH1F("fHmpt", "Relative Pt resolution (pt>4GeV/c)", 30, -60., 60.);
  99   fE = new TH1F("fE", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
 100   fEp = new TH2F("fEp", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
 101   fGoodPhi = new TH1F("fGoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
 102   fFoundPhi = new TH1F("fFoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
 103
 104   fListOfHistos->Add(fGood);
 105   fListOfHistos->Add(fFound);
 106   fListOfHistos->Add(fFake);
 107   fListOfHistos->Add(fP);
 108   fListOfHistos->Add(fL);
 109   fListOfHistos->Add(fPt);
 110   fListOfHistos->Add(fHmpt);
 111   fListOfHistos->Add(fE);
 112   fListOfHistos->Add(fEp);
 113   fListOfHistos->Add(fGoodPhi);
 114   fListOfHistos->Add(fFoundPhi);
 115 }
 116
 117
 118 void AliTPCComparisonTask::UserExec(Option_t *)
 119 {
 120   // Main loop
 121   // Called for each event
 122
 123   // MC information
 124   AliMCEvent* mcEvent = MCEvent();
 125   if (!mcEvent) {
 126     Printf("ERROR: Could not retrieve MC event");
 127     return;
 128   }
 129
 130   Int_t nt = 0;
 131
 132   TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
 133
 134   // Loop over all MC tracks and select "good" tracks
 135   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 136     AliMCParticle* track = mcEvent->GetTrack(iTracks);
 137     if (!track) {
 138       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
 139       continue;
 140     }
 141
 142     // Track selection
 143     if (track->Pt() < 0.2) continue;
 144     if (track->Pt() > 7.0) continue;
 145     if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
 146
 147     Double_t vx = track->Xv(),vy = track->Yv(),vz = track->Zv();
 148     if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
 149     if (TMath::Abs(vz) > 50.) continue;
 150
 151     // Loop over Track References
 152     Bool_t LabelTPC = kFALSE;
 153     AliTrackReference* trackRef = 0;
 154     for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
 155       trackRef = track->GetTrackReference(iTrackRef);
 156       if(trackRef) {
 157         Int_t detectorId = trackRef->DetectorId();
 158         if (detectorId == AliTrackReference::kTPC) {
 159           LabelTPC = kTRUE;
 160           break;
 161         }
 162       }
 163     }
 164     // "Good" tracks
 165     if (LabelTPC) {
 166       AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
 167       ref->SetLabel(iTracks);
 168       TParticle* particle = track->Particle();
 169       Int_t pdg = particle->GetPdgCode();
 170       ref->SetPDG(pdg);
 171       ref->SetPz(trackRef->Pz());
 172       Float_t pt = trackRef->Pt();
 173       ref->SetPt(pt);
 174       fGood->Fill(pt);
 175       Float_t phig = trackRef->Phi();
 176       ref->SetPhi(phig);
 177       fGoodPhi->Fill(phig);
 178       nt++;
 179     }
 180   } //track loop
 181
 182
 183   // ESD information
 184   AliVEvent* event = InputEvent();
 185   if (!event) {
 186     Printf("ERROR: Could not retrieve event");
 187     return;
 188   }
 189
 190   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 191
 192   Bool_t iFound;
 193   Int_t nfound = 0;
 194   Int_t nfake = 0;
 195   Int_t nlost = 0;
 196
 197   Int_t mcGoods = refs->GetEntriesFast();
 198
 199   for (Int_t k = 0; k < mcGoods; k++) {
 200     AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
 201     if (!ref) continue;
 202     Int_t mcLabel = ref->GetLabel();
 203     Float_t ptg = ref->GetPt();
 204     Float_t Phig = ref->GetPhi();
 205
 206     iFound = kFALSE;
 207     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 208       AliESDtrack* track = esd->GetTrack(iTrack);
 209       if (!track) {
 210         Printf("ERROR: Could not receive track %d", iTrack);
 211         continue;
 212       }
 213
 214       if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;
 215
 216       Int_t tpcLabel =  track->GetTPCLabel();
 217       if (mcLabel == TMath::Abs(tpcLabel)) {
 218         if (mcLabel == tpcLabel) {
 219           nfound++;
 220           fFound->Fill(ptg);
 221           fFoundPhi->Fill(Phig);
 222         }
 223         else {
 224           nfake++;
 225           fFake->Fill(ptg);
 226         }
 227         iFound = kTRUE;
 228
 229
 230         Double_t pxpypz[3]; track->GetInnerPxPyPz(pxpypz);
 231         Float_t phi = TMath::ATan2(pxpypz[1],pxpypz[0]);
 232         if (phi < 0) phi += 2*TMath::Pi();
 233         Double_t pt = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
 234         Float_t lam = TMath::ATan2(pxpypz[2],pt);
 235         Float_t pt1 = 1/pt;
 236
 237         Int_t pdg = ref->GetPDG();
 238         if (TMath::Abs(pdg) == 11 && ptg>4.) {
 239           //high momentum electrons
 240           fHmpt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);
 241         }
 242         else {
 243           fP->Fill((phi - Phig)*1000.);
 244           Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
 245           fL->Fill((lam - lamg)*1000.);
 246           fPt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);
 247         }
 248
 249         Float_t mom = pt/TMath::Cos(lam);
 250         Float_t dedx = track->GetTPCsignal();
 251         fEp->Fill(mom, dedx, 1.);
 252         if (TMath::Abs(pdg) == 211) { // pions
 253           if (mom > 0.4 && mom < 0.5) {
 254             fE->Fill(dedx,1.);
 255           }
 256         }
 257
 258         break;
 259       }
 260     }
 261     if (!iFound) {
 262       nlost++;
 263     }
 264   }
 265
 266   Printf(" Results: " );
 267   Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
 268
 269   refs->Clear();
 270
 271   // Post output data.
 272   PostData(1, fListOfHistos);
 273 }
 274
 275
 276 void AliTPCComparisonTask::Terminate(Option_t *)
 277 {
 278   // Draw result to the screen
 279   // Called once at the end of the query
 280
 281   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 282   if (!fListOfHistos) {
 283     Printf("ERROR: fListOfHistos not available");
 284     return;
 285   }
 286
 287   fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
 288   fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
 289   fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
 290   fP = dynamic_cast<TH1F*>(fListOfHistos->At(3));
 291   fL = dynamic_cast<TH1F*>(fListOfHistos->At(4));
 292   fPt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
 293   fHmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
 294   fE = dynamic_cast<TH1F*>(fListOfHistos->At(7));
 295   fEp = dynamic_cast<TH2F*>(fListOfHistos->At(8));
 296   fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
 297   fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
 298
 299   gStyle->SetOptStat(111110);
 300   gStyle->SetOptFit(1);
 301
 302   TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
 303
 304   Int_t minc = 33;
 305
 306   TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
 307   p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
 308   fP->SetFillColor(4); fP->SetXTitle("(mrad)");
 309   if (fP->GetEntries() < minc) fP->Draw(); else fP->Fit("gaus"); c1->cd();
 310
 311   TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1, 0.6); p2->Draw();
 312   p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
 313   fL->SetFillColor(4); fL->SetXTitle("(mrad)");
 314   if (fL->GetEntries() < minc) fL->Draw(); else fL->Fit("gaus"); c1->cd();
 315
 316   TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
 317   p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
 318   fPt->SetFillColor(2); fPt->SetXTitle("(%)");
 319   if (fPt->GetEntries() < minc) fPt->Draw(); else fPt->Fit("gaus"); c1->cd();
 320
 321   TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
 322   p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
 323   fHmpt->SetFillColor(6); fHmpt->SetXTitle("(%)");
 324   if (fHmpt->GetEntries() < minc) fHmpt->Draw(); else fHmpt->Fit("gaus"); c1->cd();
 325
 326   TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
 327   p5->SetFillColor(41); p5->SetFrameFillColor(10);
 328   fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
 329   TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
 330   TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
 331   hg->Divide(fFound, fGood, 1., 1., "B");
 332   hf->Divide(fFake, fGood, 1., 1., "B");
 333   hg->SetLineColor(4); hg->SetLineWidth(2);
 334   hg->SetMaximum(1.4);
 335   hg->SetYTitle("Tracking efficiency");
 336   hg->SetXTitle("Pt (GeV/c)");
 337   hg->Draw();
 338
 339   TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
 340   line1->Draw("same");
 341   TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
 342   line2->Draw("same");
 343
 344   hf->SetFillColor(1);
 345   hf->SetFillStyle(3013);
 346   hf->SetLineColor(2);
 347   hf->SetLineWidth(2);
 348   hf->Draw("histsame");
 349   TText* text = new TText(0.461176, 0.248448, "Fake tracks");
 350   text->SetTextSize(0.05);
 351   text->Draw();
 352   text = new TText(0.453919, 1.11408, "Good tracks");
 353   text->SetTextSize(0.05);
 354   text->Draw();
 355
 356   TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
 357
 358   TPad* p6 = new TPad("p6", "", 0., 0., 1., 0.5); p6->Draw();
 359   p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
 360   fE->SetFillColor(2); fE->SetFillStyle(3005);
 361   fE->SetXTitle("Arbitrary Units");
 362   if (fE->GetEntries() < minc) fE->Draw(); else fE->Fit("gaus"); c2->cd();
 363
 364   TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
 365   p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
 366   fEp->SetFillColor(2); fEp->SetFillStyle(3005);
 367   fEp->SetMarkerStyle(8);
 368   fEp->SetMarkerSize(0.4);
 369   fEp->SetXTitle("p (GeV/c)"); fEp->SetYTitle("dE/dX (Arb. Units)");
 370   fEp->Draw(); c1->cd();
 371
 372   TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
 373   c3->SetFillColor(42); c3->SetFrameFillColor(10);
 374   fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
 375   TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
 376                          90, 0., 2.*TMath::Pi());
 377   hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
 378   hgphi->SetYTitle("Tracking efficiency");
 379   hgphi->SetXTitle("Phi (rad)");
 380   hgphi->Draw();
 381
 382   TFile fc("AliTPCComparison.root", "RECREATE");
 383   c1->Write();
 384   c2->Write();
 385   c3->Write();
 386   fc.Close();
 387
 388 }