PWG1/comparison/AliITSComparisonTask.cxx

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