PWGPP/comparison/AliTRDComparisonTask.cxx

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