PWGPP/comparison/AliITSComparisonTask.cxx

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