PWG1/comparison/AliITSComparisonTask.cxx

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