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