1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /* $Id: AliTRDefficiencyMC.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Reconstruction QA //
23 // Markus Fasel <M.Fasel@gsi.de> //
25 ////////////////////////////////////////////////////////////////////////////
27 #include <TObjArray.h>
28 #include <TClonesArray.h>
33 #include <TDatabasePDG.h>
34 #include "TTreeStream.h"
38 #include "AliESDtrack.h"
39 #include "AliMathBase.h"
40 #include "AliTrackReference.h"
41 #include "AliAnalysisManager.h"
43 #include "AliTRDcluster.h"
44 #include "AliTRDseedV1.h"
45 #include "AliTRDtrackV1.h"
46 #include "AliTRDpidUtil.h"
47 #include "Cal/AliTRDCalPID.h"
48 #include "info/AliTRDtrackInfo.h"
49 #include "AliTRDinfoGen.h"
50 #include "AliTRDefficiencyMC.h"
52 ClassImp(AliTRDefficiencyMC)
53 Float_t AliTRDefficiencyMC::fgPCut = 0.2; //[GeV/c]
54 Float_t AliTRDefficiencyMC::fgPhiCut = 50.; //[deg]
55 Float_t AliTRDefficiencyMC::fgThtCut = 50.; //[deg]
56 //_____________________________________________________________________________
57 AliTRDefficiencyMC::AliTRDefficiencyMC()
61 // Default constructor
65 AliTRDefficiencyMC::AliTRDefficiencyMC(char* name)
66 :AliTRDrecoTask(name, "Combined Tracking Efficiency")
69 // Default constructor
73 //_____________________________________________________________________________
74 void AliTRDefficiencyMC::UserCreateOutputObjects(){
76 // Create output objects
79 fContainer = Histos();
80 PostData(1, fContainer);
83 //_____________________________________________________________________________
84 void AliTRDefficiencyMC::UserExec(Option_t *){
88 // Loop over TrackInfos
89 // 1st: check if there is a trackTRD
90 // 2nd: put conditions on the track:
91 // - check if we did not register it before
92 // - check if we have Track References for the track
93 // 3rd: Register track:
94 // - accepted if both conditions are fulfilled
95 // - contamination if at least one condition is not fulfilled
96 // 4th: check Monte-Carlo Track wheter findable or not if there is no TRD track in this track info
97 // 5th: register MC track as rejected if findable and not jet registered
98 // Match the registers accepted and rejected and clean register rejected
99 // Fill the histograms
101 const Int_t kArraySize = 10000; // Store indices of track references in an array
102 Int_t indexAccept[kArraySize],
103 indexReject[kArraySize],
104 indexContam[kArraySize];
105 memset(indexAccept, 0, sizeof(Int_t) * kArraySize);
106 memset(indexReject, 0, sizeof(Int_t) * kArraySize);
107 memset(indexContam, 0, sizeof(Int_t) * kArraySize);
113 Bool_t isContamination = kFALSE;
115 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
117 Int_t nTrackInfos = fTracks->GetEntriesFast();
118 AliDebug(2, Form(" CANDIDATE TRACKS[%d]", nTrackInfos));
120 AliTRDtrackV1 *trackTRD(NULL);
121 AliTRDtrackInfo *trkInf(NULL);
122 for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
123 trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
124 if(!trkInf) continue;
126 if(trkInf->GetTrack() || trkInf->GetNumberOfClustersRefit()){
127 isContamination = (IsRegistered(trkInf,indexAccept,naccept)>=0);
128 if(!trkInf->GetNTrackRefs()){
129 // We reject the track since the Monte Carlo Information is missing
130 AliDebug(2, Form("MC(Track Reference) missing @ label[%d]", trkInf->GetLabel()));
131 isContamination = kTRUE;
133 if(trackTRD && DebugLevel()>5) FillStreamTrackWOMC(trkInf);
136 // reject kink (we count these only once)
137 if(trkInf->GetKinkIndex()){
138 AliDebug(4, Form(" track @ idx[%d] MC[%d] is kink.", itinf, trkInf->GetLabel()));
142 // Register track as contamination
143 AliDebug(4, Form(" track @ idx[%d] MC[%d] is contamination.", itinf, trkInf->GetLabel()));
144 indexContam[ncontam++]=itinf;
148 AliDebug(4, Form(" track @ idx[%d] is ACCEPTED.", itinf));
149 // Register track as accepted
150 indexAccept[naccept++] = itinf;
153 if((code=IsFindableNot(trkInf))){
154 AliDebug(4, Form(" track @ idx[%d] MC[%d] not findable [%d].", itinf, trkInf->GetLabel(), code));
157 // register track as rejected if not already registered there
159 // These track infos are not!!! registered as contamination
160 if(IsRegistered(trkInf, indexReject, nreject)<0){
161 AliDebug(4, Form(" track @ idx[%d] MC[%d] is missed.", itinf, trkInf->GetLabel()));
162 indexReject[nreject++] = itinf;
167 AliDebug(2, Form("TRACKS STATISTICS naccept[%d] ncontam[%d] nkink[%d] nmissed[%d] nfindnt[%d] ALL[%d] LOST[%d]", naccept, ncontam, nkink, nreject, nfindnt, naccept+ncontam+nkink, nTrackInfos-(naccept+nreject+ncontam+nkink+nfindnt)));
169 // we have to check if the rejected tracks are registered as found
170 // a possible source for this:
171 // first the track is not found by the barrel tracking but it is later found
172 // by the stand alone tracking, then two track info objects with the same
173 // label would be created
175 // these tracks are not! registered as contamination
176 Int_t tmprejected[kArraySize]; Int_t nrej = nreject;
177 memcpy(tmprejected, indexReject, sizeof(Int_t) * nreject);
179 for(Int_t irej = 0; irej < nrej; irej++){
180 trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej]));
182 if((idx = IsRegistered(trkInf,indexAccept,naccept))<0){
183 indexReject[nreject++] = tmprejected[irej];
185 //printf("tracks @ accept[%d] missed[%d] are the same.\n", indexAccept[idx], tmprejected[irej]);
190 FillHistograms(naccept, &indexAccept[0], kAccept);
191 FillHistograms(nreject, &indexReject[0], kMiss);
192 FillHistograms(ncontam, &indexContam[0], kFake);
194 Int_t nall(naccept + nreject);
195 AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] Fake[%3d | %5.2f%%]",
196 (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(),
199 (nall ? 1.E2*Float_t(naccept)/Float_t(nall) : 0.),
201 (nall ? 1.E2*Float_t(ncontam)/Float_t(nall) : 0.)));
205 //_____________________________________________________________________________
206 Bool_t AliTRDefficiencyMC::PostProcess()
211 // Change the histogram style
212 // For species histograms apply the colors connected with the given particle species
217 //_____________________________________________________________________________
218 Bool_t AliTRDefficiencyMC::GetRefFigure(Int_t ifig){
220 // Plot the histograms
222 if(ifig >= fNRefFigures) return kFALSE;
223 if(!gPad) return kFALSE;
224 gPad->SetLogx(kTRUE);
226 (dynamic_cast<TH1 *>(fContainer->At(ifig)))->Draw("e1");
230 TLegend *leg=new TLegend(.65, .12, .85, .3);
231 leg->SetHeader("Charge");
232 leg->SetBorderSize(1);leg->SetFillColor(kWhite);
235 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram));
236 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
237 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+1));
238 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
242 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+2));
243 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
244 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+3));
245 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
249 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+4));
250 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
251 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+5));
252 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
256 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+6));
257 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
258 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+7));
259 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
263 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+8));
264 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
265 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+9));
266 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
270 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+10));
271 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
272 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+11));
273 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
277 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+12));
278 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
279 h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+13));
280 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
287 //_____________________________________________________________________________
288 TObjArray *AliTRDefficiencyMC::Histos(){
290 // Create the histograms
293 if(fContainer) return fContainer;
294 const Int_t nbins = AliTRDCalPID::kNMom;
295 Float_t xbins[nbins+1] = {fgPCut, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
296 const Int_t marker[2][AliPID::kSPECIES+1] = {
297 {20, 21, 22, 23, 29, 2},
298 {24, 25, 26, 27, 30, 5}
301 fContainer = new TObjArray();fContainer->Expand(14);
304 fContainer->AddAt(h=new TProfile("hEff", "Tracking Efficiency ALL", nbins, xbins), kEfficiencyHistogram);
305 h->SetMarkerStyle(22);
306 h->SetMarkerColor(kBlue);
307 h->GetXaxis()->SetTitle("p [GeV/c]");
308 h->GetXaxis()->SetMoreLogLabels();
309 h->GetYaxis()->SetTitle("Efficiency");
310 h->GetYaxis()->SetRangeUser(0.2, 1.1);
311 fContainer->AddAt(h=new TProfile("hFake", "Fake Tracks", nbins, xbins), kContaminationHistogram);
312 h->SetMarkerStyle(22);
313 h->SetMarkerColor(kBlue);
314 h->GetXaxis()->SetTitle("p [GeV/c]");
315 h->GetXaxis()->SetMoreLogLabels();
316 h->GetYaxis()->SetTitle("Contamination");
318 Char_t sign[]={'+', '-'};
319 for(Int_t isign = 0; isign < 2; isign++){
320 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
321 fContainer->AddAt(h=new TProfile(
322 Form("hEff_%s%c", AliPID::ParticleShortName(ispec), sign[isign]),
323 Form("Tracking Efficiency for %s", AliPID::ParticleName(ispec)), nbins, xbins),
324 kEfficiencySpeciesHistogram+ispec*2+isign);
325 h->SetMarkerStyle(marker[isign][ispec]);
326 h->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
327 h->SetMarkerColor(kBlack);
328 h->GetXaxis()->SetTitle("p [GeV/c]");
329 h->GetXaxis()->SetMoreLogLabels();
330 h->GetYaxis()->SetTitle("Efficiency");
331 h->GetYaxis()->SetRangeUser(0.2, 1.1);
334 fContainer->AddAt(h=new TProfile(Form("hEff_PID%c", sign[isign]), "Tracking Efficiency no PID", nbins, xbins), kEfficiencySpeciesHistogram+AliPID::kSPECIES*2+isign);
335 h->SetMarkerStyle(marker[isign][AliPID::kSPECIES]);
336 h->SetMarkerColor(kBlack);h->SetLineColor(kBlack);
337 h->GetXaxis()->SetTitle("p [GeV/c]");
338 h->GetXaxis()->SetMoreLogLabels();
339 h->GetYaxis()->SetTitle("Efficiency");
340 h->GetYaxis()->SetRangeUser(0.2, 1.1);
345 //_____________________________________________________________________________
346 Int_t AliTRDefficiencyMC::IsFindableNot(AliTRDtrackInfo * const trkInf){
348 // Apply Cuts on the Monte Carlo track references
349 // return whether track is findable or not
353 const Float_t chmbHght = AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
354 const Float_t eps(1.E-3);
355 Int_t ntr(trkInf->GetNTrackRefs());
357 AliDebug(10, Form(" CANDIDATE TrackRefs[%d]", ntr));
358 // Check if track is findable
359 Double_t mom(0.), phi(0.), tht(0.);
360 Float_t xmin = 10000.0, xmax = 0.0;
361 Float_t ymin = 0.0, ymax = 0.0;
362 Float_t zmin = 0.0, zmax = 0.0;
363 Float_t lastx = 0.0, x = 0.0;
364 Int_t nLayers(0), ntrTRD(0);
366 AliTrackReference *trackRef(NULL);
367 for(Int_t itr(0); itr<ntr; itr++){
368 if(!(trackRef = trkInf->GetTrackRef(itr))) continue;
369 x = trackRef->LocalX();
370 // Be Sure that we are inside TRD
371 if(x < AliTRDinfoGen::GetEndTPC() || x > AliTRDinfoGen::GetEndTRD()) continue;
372 sector[ntrTRD] = Int_t(trackRef->Alpha()/AliTRDgeometry::GetAlpha());
373 AliDebug(10, Form(" [%2d] x[%7.2f] y[%7.2f] z[%7.2f] Sec[%2d]", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z(), sector[ntrTRD]));
375 xmin = trackRef->LocalX();
376 ymin = trackRef->LocalY();
377 zmin = trackRef->Z();
380 xmax = trackRef->LocalX();
381 ymax = trackRef->LocalY();
382 zmax = trackRef->Z();
385 Float_t dist = TMath::Abs(x - lastx);
386 if(TMath::Abs(dist - chmbHght) < eps && sector[ntrTRD]==sector[0]){
387 AliDebug(10, Form(" dx = %7.2f", dist));
392 ntrTRD++; if(ntrTRD>=20) break;
394 Double_t dx(xmax - xmin);
395 if(TMath::Abs(dx)<eps) return kNoChmb;
397 phi = (ymax -ymin)/dx;
398 tht = (zmax -zmin)/dx;
399 phi=TMath::ATan(phi)*TMath::RadToDeg();
400 tht=TMath::ATan(tht)*TMath::RadToDeg();
401 Bool_t primary = trkInf->IsPrimary();
402 const AliTRDtrackInfo::AliESDinfo *esd(trkInf->GetESDinfo());
403 AliDebug(10, Form(" p=%6.3f[GeV/c] phi=%6.2f[deg] theta=%6.2f[deg] nLy[%d]",
404 mom, phi, tht, nLayers));
406 (*DebugStream()) << "IsFindable"
412 << "Primary=" << primary
417 if(!nLayers) return kNoChmb;
418 if(xmax < xmin) return kCurved;
419 if(mom < fgPCut) return kPCut;
422 if(TMath::Abs(phi) > fgPhiCut) return kPhiCut;
423 if(TMath::Abs(tht) > fgThtCut) return kThtCut;
426 if(!esd)return kLayer;
427 if(!(esd->GetStatus() & AliESDtrack::kTPCout)) return kLayer;
430 //if(!trkInf->IsPrimary()) {failCode=kPrimary; return kFALSE;}
435 //_____________________________________________________________________________
436 void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, ETRDefficiencyMCstatus mode){
438 // Fill Histograms in three different modes:
439 // 1st tracks which are found and accepted
440 // 2nd tracks which are not found and not already counted
441 // 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
444 TDatabasePDG *dbPDG(TDatabasePDG::Instance());
445 Double_t trkmom(0.); // the track momentum
446 Int_t trkpdg(-1); // particle PDG code
447 AliTRDtrackInfo *trkInf(NULL);
448 for(Int_t itk = 0; itk < nTracks; itk++){
449 trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk]));
450 if(trkInf->GetNTrackRefs()){
451 // use Monte-Carlo Information for Momentum and PID
452 trkmom = trkInf->GetTrackRef(0)->P();
453 trkpdg = trkInf->GetPDG();
456 trkmom = trkInf->GetTrack()->P();
459 const Char_t *cmode(NULL);
462 (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 1);
463 (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
467 (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 0);
468 (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
472 (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 1);
476 AliDebug(3, Form(" track[%d] MC[%d] Mode[%s]", indices[itk], trkInf->GetLabel(), cmode));
478 // Fill species histogram
479 Int_t idxSpec = AliTRDpidUtil::Pdg2Pid(TMath::Abs(trkpdg));
480 Int_t sign = dbPDG->GetParticle(trkpdg)->Charge() > 0. ? 1 : 0;
481 //printf("[%d]%s pdg[%d] sign[%d]\n", idxSpec, AliPID::ParticleName(idxSpec), trkpdg, sign);
482 if(idxSpec < 0) idxSpec = AliPID::kSPECIES;
483 (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + idxSpec*2+sign)))->Fill(trkmom, mode==kAccept?1:0);
487 //_____________________________________________________________________________
488 void AliTRDefficiencyMC::FillStreamTrackWOMC(AliTRDtrackInfo * const trkInf){
491 // 1. The event number
492 // 2. The track label
493 // 3. The TRD track label
494 // 4. The frequency of the TRD Label
495 // 5. Momentum from TPC (NO STAND ALONE TRACK)
498 // 8. Monte Carlo PID
499 // 9. We check the Labels of the TRD track according to them we search the maching Monte-Carlo track.
500 // From the matching Monte-Carlo track we store trackRefs, phi and momentum
501 // 10. We may also want to keep the kink index
502 Double_t mom = trkInf->GetESDinfo()->GetOuterParam()->P();
503 Int_t event = (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry();
504 Int_t label = trkInf->GetLabel();
505 Int_t kinkIndex = trkInf->GetKinkIndex();
506 Int_t pdg = trkInf->GetPDG();
507 Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
508 Int_t labelsTRD[180]; // Container for the cluster labels
509 Int_t sortlabels[360]; // Cluster Labels sorted according their occurancy
510 AliTRDseedV1 *tracklet(NULL);
511 AliTRDcluster *c(NULL);
513 AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
514 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
515 tracklet = trackTRD->GetTracklet(il);
516 if(!tracklet) continue;
517 tracklet->ResetClusterIter();
519 while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
521 // Determine Label and Frequency
522 AliMathBase::Freq(nclusters, const_cast<const Int_t *>(&labelsTRD[0]), &sortlabels[0], kTRUE);
523 Int_t labelTRD = sortlabels[0];
524 Int_t freqTRD = sortlabels[1];
525 // find the track info object matching to the TRD track
526 AliTRDtrackInfo *realtrack = 0;
527 TObjArrayIter rtiter(fTracks);
528 while((realtrack = (AliTRDtrackInfo *)rtiter())){
529 if(realtrack->GetLabel() != labelTRD) continue;
532 TClonesArray trackRefs("AliTrackReference");
535 Double_t realPhi = 0.;
537 // pack the track references into the trackRefsContainer
538 for(Int_t iref = 0; iref < realtrack->GetNTrackRefs(); iref++){
539 new(trackRefs[iref])AliTrackReference(*(realtrack->GetTrackRef(iref)));
541 realPdg = realtrack->GetPDG();
542 if(realtrack->GetNTrackRefs()){
543 realP = realtrack->GetTrackRef(0)->P();
544 realPhi = realtrack->GetTrackRef(0)->Phi();
547 (*DebugStream()) << "EffMCfake"
550 << "labelTRD=" << labelTRD
551 << "FreqTRDlabel=" << freqTRD
553 << "phiTPC=" << phiTPC
554 << "trackTRD=" << trackTRD
556 << "TrackRefs=" << &trackRefs
557 << "RealPDG=" << realPdg
559 << "RealPhi" << realPhi
560 << "KinkIndex=" << kinkIndex
564 //_____________________________________________________________________________
565 Int_t AliTRDefficiencyMC::IsRegistered(AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
567 // Checks if track is registered in a given mode
570 Int_t label(trkInf->GetLabel());
571 for(Int_t il(nTracks); il--;){
572 if((dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))->GetLabel() == label) return il;