PWGPP/comparison/AliTOFComparisonTask.cxx

   1 //-------------------------------------------------------------------------
   2 //
   3 // This is the PROOF-enabled version of TOF/AliTOFComparison.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
  28 #include "AliTOFComparisonTask.h"
  29
  30 ClassImp(AliTOFComparisonTask)
  31
  32 extern TStyle *gStyle;
  33
  34 AliTOFComparisonTask::AliTOFComparisonTask()
  35   : AliAnalysisTaskSE("AliTOFComparisonTask"),
  36     fListOfHistos(0),
  37     fGood(0),
  38     fFound(0),
  39     fFake(0),
  40     fGoodPhi(0),
  41     fFoundPhi(0),
  42     fGoodl(0),
  43     fFakel(0),
  44     fFoundl(0)
  45 {
  46   // Default constructor
  47   AliInfo("Defaulut Constructor AliTOFComparisonTask");
  48   // Define input and output slots here
  49   // Input slot #0 works with a TChain
  50   DefineInput(0, TChain::Class());
  51   // Output slot #1 TList
  52   DefineOutput(1, TList::Class());
  53 }
  54
  55 AliTOFComparisonTask::AliTOFComparisonTask(const char* name)
  56   : AliAnalysisTaskSE(name),
  57     fListOfHistos(0),
  58     fGood(0),
  59     fFound(0),
  60     fFake(0),
  61     fGoodPhi(0),
  62     fFoundPhi(0),
  63     fGoodl(0),
  64     fFakel(0),
  65     fFoundl(0)
  66 {
  67   // Constructor
  68   AliInfo("Constructor AliTOFComparisonTask");
  69   // Define input and output slots here
  70   // Input slot #0 works with a TChain
  71   DefineInput(0, TChain::Class());
  72   // Output slot #1 TList
  73   DefineOutput(1, TList::Class());
  74 }
  75
  76
  77
  78 void AliTOFComparisonTask::UserCreateOutputObjects()
  79 {
  80   // Create histograms
  81   // Called once
  82   AliInfo("AliTOFComparisonTask::UserCreateOutputObjects");
  83   // Create output container
  84   fListOfHistos = new TList();
  85
  86   fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
  87   fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
  88   fFake = new TH1F("fFake",  "Pt for fake tracks", 34, 0.2, 7.0);
  89   fGoodPhi = new TH1F("fGoodPhi", "Phi for Good tracks", 90, 0., 2.*TMath::Pi());
  90   fFoundPhi = new TH1F("fFoundPhi", "Phi for Found tracks", 90, 0., 2.*TMath::Pi());
  91   fGoodl = new TH1F("fGoodl", "Good tracks", 30, -1., 1.);
  92   fFakel = new TH1F("fFakel", "Mismatched tracks", 30, -1., 1.);
  93   fFoundl = new TH1F("fFoundl", "Matched tracks", 30, -1., 1.);
  94
  95   fListOfHistos->Add(fGood);
  96   fListOfHistos->Add(fFound);
  97   fListOfHistos->Add(fFake);
  98   fListOfHistos->Add(fGoodPhi);
  99   fListOfHistos->Add(fFoundPhi);
 100   fListOfHistos->Add(fGoodl);
 101   fListOfHistos->Add(fFakel);
 102   fListOfHistos->Add(fFoundl);
 103 }
 104
 105
 106 void AliTOFComparisonTask::UserExec(Option_t *)
 107 {
 108   // Main loop
 109   // Called for each event
 110
 111   // MC information
 112   AliMCEvent* mcEvent = MCEvent();
 113   if (!mcEvent) {
 114     Printf("ERROR: Could not retrieve MC event");
 115     return;
 116   }
 117
 118   Int_t nt = 0;
 119
 120   TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
 121
 122   // Loop over all MC tracks and select "good" tracks
 123   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 124     AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
 125     if (!track) {
 126       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
 127       continue;
 128     }
 129
 130     // Track selection
 131     if (track->Pt() < 0.2) continue;
 132     if (track->Pt() > 7.0) continue;
 133     if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
 134
 135     Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
 136     if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
 137     if (TMath::Abs(vz) > 50.) continue;
 138
 139     if (TMath::Abs(vx) > 0.1) continue;
 140     if (TMath::Abs(vy) > 0.1) continue;
 141
 142
 143     // Loop over Track References
 144     Bool_t labelTPC = kFALSE;
 145     Bool_t labelTOF = kFALSE;
 146     AliTrackReference* trackRef = 0;
 147     UInt_t iTSLayerMap = 0;
 148     for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
 149       trackRef = track->GetTrackReference(iTrackRef);
 150       if(trackRef) {
 151         Int_t detectorId = trackRef->DetectorId();
 152         if (detectorId == AliTrackReference::kITS) {
 153           Float_t Radius = trackRef->R();
 154           if (Radius > 2.5 && Radius < 5.5)
 155             iTSLayerMap |= 1;
 156           else if (Radius > 5.5 && Radius < 8.5)
 157             iTSLayerMap |= 1 << 1;
 158           else if (Radius > 13.5 && Radius < 16.5)
 159             iTSLayerMap |= 1 << 2;
 160           else if (Radius > 22. && Radius < 26.)
 161             iTSLayerMap |= 1 << 3;
 162           else if (Radius > 36. && Radius < 41.)
 163             iTSLayerMap |= 1 << 4;
 164           else if (Radius > 42. && Radius < 46.)
 165             iTSLayerMap |= 1 << 5;
 166           else {
 167             Printf("Wrong radius %f ", Radius);
 168             return;
 169           }
 170
 171         }
 172         if (detectorId == AliTrackReference::kTPC) {
 173           labelTPC = kTRUE;
 174         }
 175         if (detectorId == AliTrackReference::kTOF) {
 176           labelTOF = kTRUE;
 177         }
 178       }
 179     }    // track references loop
 180
 181     // Skip tracks that passed not all ITS layers
 182     if (iTSLayerMap != 0x3F) continue;
 183
 184     // "Good" tracks
 185     if (labelTPC && labelTOF) {
 186       AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
 187       ref->SetLabel(iTracks);
 188       TParticle* particle = track->Particle();
 189       Int_t pdg = particle->GetPdgCode();
 190       ref->SetPDG(pdg);
 191       ref->SetPz(track->Pz());
 192       Float_t Pt = track->Pt();
 193       ref->SetPt(Pt);
 194       fGood->Fill(Pt);
 195       Float_t phig = track->Phi();
 196       ref->SetPhi(phig);
 197       fGoodPhi->Fill(phig);
 198       Double_t tgl = track->Pz()/Pt;
 199       fGoodl->Fill(tgl);
 200       nt++;
 201     }
 202   }    // track loop
 203
 204   // ESD information
 205   AliVEvent* event = InputEvent();
 206   if (!event) {
 207     Printf("ERROR: Could not retrieve event");
 208     return;
 209   }
 210
 211   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 212
 213   Bool_t iFound;
 214   Int_t nfound = 0;
 215   Int_t nfake = 0;
 216   Int_t nlost = 0;
 217
 218   Int_t mcGoods = refs->GetEntriesFast();
 219
 220   // Loop over all "good" MC tracks
 221   for (Int_t k = 0; k < mcGoods; k++) {
 222     AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
 223     if (!ref) continue;
 224     Int_t MCLabel = ref->GetLabel();
 225     Float_t ptg = ref->GetPt();
 226     Float_t Phig = ref->GetPhi();
 227     iFound = kFALSE;
 228     Int_t iTrack;
 229     AliESDtrack* track = 0;
 230     for (iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 231       track = esd->GetTrack(iTrack);
 232       if (!track) {
 233         Printf("ERROR: Could not receive track %d", iTrack);
 234         continue;
 235       }
 236
 237       Int_t Label =  track->GetLabel();
 238       if (MCLabel != TMath::Abs(Label)) continue;
 239
 240       if (track->GetTOFsignal() <= 0.) continue;
 241
 242       Double_t tgl = ref->GetPz()/ref->GetPt();
 243       Int_t TOFLabel[3];
 244       track->GetTOFLabel(TOFLabel);
 245       if (TOFLabel[0] != MCLabel)
 246         if (TOFLabel[1] != MCLabel)
 247           if (TOFLabel[2] != MCLabel) {
 248             fFake->Fill(ptg);
 249             fFakel->Fill(tgl);
 250             nfake++;
 251             break;
 252           }
 253       fFound->Fill(ptg);
 254       fFoundl->Fill(tgl);
 255       fFoundPhi->Fill(Phig);
 256       nfound++;
 257       break;
 258     }
 259     if(iTrack == esd->GetNumberOfTracks()) {
 260       Printf("Not matched: MCLabel %d ", MCLabel);
 261       nlost++;
 262       if (track) {
 263         Printf("  kITSout %d kTPCout %d kTRDout %d kTIME %d",
 264                (track->GetStatus()&AliESDtrack::kITSout),
 265                (track->GetStatus()&AliESDtrack::kTPCout),
 266                (track->GetStatus()&AliESDtrack::kTRDout),
 267                (track->GetStatus()&AliESDtrack::kTIME));
 268       }
 269       else
 270         Printf(" No ESD track !");
 271     }
 272   }
 273
 274   Printf(" Results: " );
 275   Printf(" Good %d Found %d Fake %d Lost %d", nt, nfound, nfake, nlost);
 276
 277   refs->Clear();
 278
 279   // Post output data.
 280   PostData(1, fListOfHistos);
 281
 282 }
 283
 284
 285 void AliTOFComparisonTask::Terminate(Option_t *)
 286 {
 287   // Draw result to the screen
 288   // Called once at the end of the query
 289   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 290   if (!fListOfHistos) {
 291     Printf("ERROR: fListOfHistos not available");
 292     return;
 293   }
 294
 295   fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
 296   fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
 297   fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
 298   fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3));
 299   fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(4));
 300   fGoodl = dynamic_cast<TH1F*>(fListOfHistos->At(5));
 301   fFakel = dynamic_cast<TH1F*>(fListOfHistos->At(6));
 302   fFoundl = dynamic_cast<TH1F*>(fListOfHistos->At(7));
 303
 304   gStyle->SetOptStat(111110);
 305   gStyle->SetOptFit(1);
 306
 307   TCanvas* c1 = new TCanvas("c1","",0,0,600,900);
 308
 309   TPad* p1 = new TPad("p1", "", 0., 0.5, 1., 1.); p1->Draw();
 310   p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
 311
 312   fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
 313   TH1F* hgp = new TH1F("hgp", " ", 34, 0.2, 7.0);
 314   TH1F* hfp = new TH1F("hfp", "Probability of mismatching", 34, 0.2, 7.0);
 315   hgp->Divide(fFound, fGood, 1., 1., "b");
 316   hfp->Divide(fFake, fGood, 1., 1., "b");
 317   hgp->SetLineColor(4); hgp->SetLineWidth(2);
 318   hgp->SetMaximum(1.4);
 319   hgp->SetYTitle("Tracking efficiency");
 320   hgp->SetXTitle("Pt (GeV/c)");
 321   hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
 322   hfp->SetLineColor(2);
 323
 324   hgp->Draw();
 325   hfp->Draw("histsame");
 326   TLine *line1 = new TLine(0.2,1.0,7.0,1.0); line1->SetLineStyle(4);
 327   line1->Draw("same");
 328   TLine *line2 = new TLine(0.2,0.9,7.0,0.9); line2->SetLineStyle(4);
 329   line2->Draw("same");
 330   c1->cd();
 331
 332   TPad* p2 = new TPad("p2", "", 0., 0., 1., 0.5); p2->Draw();
 333   p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
 334
 335   fFoundl->Sumw2(); fGoodl->Sumw2(); fFakel->Sumw2();
 336   TH1F* hgl = new TH1F("hgl", "", 30, -1., 1.);
 337   TH1F* hfl = new TH1F("hfl", "Probability of mismatching", 30, -1., 1.);
 338   hgl->Divide(fFoundl, fGoodl, 1., 1., "b");
 339   hfl->Divide(fFakel, fGoodl, 1., 1., "b");
 340   hgl->SetLineColor(4); hgl->SetLineWidth(2);
 341   hgl->SetMaximum(1.4);
 342   hgl->SetYTitle("Tracking efficiency");
 343   hgl->SetXTitle("Tan(lambda)");
 344   hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
 345   hfl->SetLineColor(2);
 346
 347   hgl->Draw();
 348   hfl->Draw("histsame");
 349   TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
 350   line3->Draw("same");
 351   TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
 352   line4->Draw("same");
 353
 354   c1->Update();
 355
 356   TCanvas* c2 = new TCanvas("c2", "", 10, 10, 510, 510);
 357   c2->SetFillColor(42); c2->SetFrameFillColor(10);
 358   fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
 359   TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
 360                          90, 0., 2.*TMath::Pi());
 361   hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
 362   hgphi->SetYTitle("Tracking efficiency");
 363   hgphi->SetXTitle("Phi (rad)");
 364   hgphi->Draw();
 365
 366   TFile fc("AliTOFComparison.root","RECREATE");
 367   c1->Write();
 368   c2->Write();
 369   fc.Close();
 370 }