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