]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDefficiencyMC.cxx
General macro for QA checks
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDefficiencyMC.cxx
CommitLineData
1ee39b3a 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
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**************************************************************************/
15
16/* $Id: AliTRDefficiencyMC.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Reconstruction QA //
21// //
22// Authors: //
23// Markus Fasel <M.Fasel@gsi.de> //
24// //
25////////////////////////////////////////////////////////////////////////////
26
27#include <TObjArray.h>
28#include <TClonesArray.h>
770382d9 29#include <TPad.h>
4226db3e 30#include <TLegend.h>
1ee39b3a 31#include <TProfile.h>
32#include <TMath.h>
770382d9 33#include <TDatabasePDG.h>
1ee39b3a 34#include "TTreeStream.h"
35
36#include "AliMagF.h"
37#include "AliPID.h"
770382d9 38#include "AliESDtrack.h"
1ee39b3a 39#include "AliMathBase.h"
40#include "AliTrackReference.h"
41#include "AliAnalysisManager.h"
42
43#include "AliTRDcluster.h"
44#include "AliTRDseedV1.h"
45#include "AliTRDtrackV1.h"
068e2c00 46#include "AliTRDpidUtil.h"
1ee39b3a 47#include "Cal/AliTRDCalPID.h"
48#include "info/AliTRDtrackInfo.h"
49#include "AliTRDinfoGen.h"
50#include "AliTRDefficiencyMC.h"
51
52ClassImp(AliTRDefficiencyMC)
770382d9 53Float_t AliTRDefficiencyMC::fgPCut = 0.2; //[GeV/c]
54Float_t AliTRDefficiencyMC::fgPhiCut = 50.; //[deg]
55Float_t AliTRDefficiencyMC::fgThtCut = 50.; //[deg]
1ee39b3a 56//_____________________________________________________________________________
57AliTRDefficiencyMC::AliTRDefficiencyMC()
f8f46e4d 58 :AliTRDrecoTask()
1ee39b3a 59{
60 //
61 // Default constructor
62 //
63}
64
f8f46e4d 65AliTRDefficiencyMC::AliTRDefficiencyMC(char* name)
66 :AliTRDrecoTask(name, "Combined Tracking Efficiency")
67{
68 //
69 // Default constructor
70 //
71}
1ee39b3a 72
73//_____________________________________________________________________________
f8f46e4d 74void AliTRDefficiencyMC::UserCreateOutputObjects(){
1ee39b3a 75 //
76 // Create output objects
77 //
78
1ee39b3a 79 fContainer = Histos();
068e2c00 80 PostData(1, fContainer);
1ee39b3a 81}
82
83//_____________________________________________________________________________
f8f46e4d 84void AliTRDefficiencyMC::UserExec(Option_t *){
1ee39b3a 85 //
86 // Execute the task:
87 //
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
100 //
101 const Int_t kArraySize = 10000; // Store indices of track references in an array
770382d9 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);
108 Int_t naccept(0),
109 nreject(0),
110 nfindnt(0),
111 nkink(0),
112 ncontam(0);
1ee39b3a 113 Bool_t isContamination = kFALSE;
114
5935a6da 115 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
116 if(!fTracks) return;
1ee39b3a 117 Int_t nTrackInfos = fTracks->GetEntriesFast();
770382d9 118 AliDebug(2, Form(" CANDIDATE TRACKS[%d]", nTrackInfos));
119
7fe4e88b 120 //AliTRDtrackV1 *trackTRD(NULL);
770382d9 121 AliTRDtrackInfo *trkInf(NULL);
1ee39b3a 122 for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
123 trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
124 if(!trkInf) continue;
770382d9 125
1ee39b3a 126 if(trkInf->GetTrack() || trkInf->GetNumberOfClustersRefit()){
770382d9 127 isContamination = (IsRegistered(trkInf,indexAccept,naccept)>=0);
1ee39b3a 128 if(!trkInf->GetNTrackRefs()){
129 // We reject the track since the Monte Carlo Information is missing
770382d9 130 AliDebug(2, Form("MC(Track Reference) missing @ label[%d]", trkInf->GetLabel()));
1ee39b3a 131 isContamination = kTRUE;
132 // Debugging
7fe4e88b 133 if(trkInf && DebugLevel()>5) FillStreamTrackWOMC(trkInf);
770382d9 134 }
1ee39b3a 135 if(isContamination){
136 // reject kink (we count these only once)
770382d9 137 if(trkInf->GetKinkIndex()){
138 AliDebug(4, Form(" track @ idx[%d] MC[%d] is kink.", itinf, trkInf->GetLabel()));
139 nkink++;
140 continue;
141 }
1ee39b3a 142 // Register track as contamination
770382d9 143 AliDebug(4, Form(" track @ idx[%d] MC[%d] is contamination.", itinf, trkInf->GetLabel()));
144 indexContam[ncontam++]=itinf;
1ee39b3a 145 continue;
146 }
147 // Accept track
770382d9 148 AliDebug(4, Form(" track @ idx[%d] is ACCEPTED.", itinf));
1ee39b3a 149 // Register track as accepted
770382d9 150 indexAccept[naccept++] = itinf;
1ee39b3a 151 }else{
770382d9 152 Int_t code(0);
153 if((code=IsFindableNot(trkInf))){
154 AliDebug(4, Form(" track @ idx[%d] MC[%d] not findable [%d].", itinf, trkInf->GetLabel(), code));
155 nfindnt++;
156 } else {
1ee39b3a 157 // register track as rejected if not already registered there
158 // Attention:
159 // These track infos are not!!! registered as contamination
770382d9 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;
163 }
1ee39b3a 164 }
165 }
166 }
770382d9 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)));
168
1ee39b3a 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
174 // Attention:
175 // these tracks are not! registered as contamination
770382d9 176 Int_t tmprejected[kArraySize]; Int_t nrej = nreject;
177 memcpy(tmprejected, indexReject, sizeof(Int_t) * nreject);
178 nreject = 0;
1ee39b3a 179 for(Int_t irej = 0; irej < nrej; irej++){
7fe4e88b 180 if(!(trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej])))) continue;
770382d9 181 Int_t idx(-1);
182 if((idx = IsRegistered(trkInf,indexAccept,naccept))<0){
183 indexReject[nreject++] = tmprejected[irej];
184 }else{
185 //printf("tracks @ accept[%d] missed[%d] are the same.\n", indexAccept[idx], tmprejected[irej]);
186 }
1ee39b3a 187 }
770382d9 188
1ee39b3a 189 // Fill Histograms
770382d9 190 FillHistograms(naccept, &indexAccept[0], kAccept);
191 FillHistograms(nreject, &indexReject[0], kMiss);
192 FillHistograms(ncontam, &indexContam[0], kFake);
193
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(),
197 nall,
198 naccept,
199 (nall ? 1.E2*Float_t(naccept)/Float_t(nall) : 0.),
200 ncontam,
201 (nall ? 1.E2*Float_t(ncontam)/Float_t(nall) : 0.)));
1ee39b3a 202}
203
204
205//_____________________________________________________________________________
206Bool_t AliTRDefficiencyMC::PostProcess()
207{
208 //
209 // Post Process
210 //
211 // Change the histogram style
212 // For species histograms apply the colors connected with the given particle species
770382d9 213 fNRefFigures = 8;
1ee39b3a 214 return kTRUE;
215}
216
217//_____________________________________________________________________________
218Bool_t AliTRDefficiencyMC::GetRefFigure(Int_t ifig){
219 //
220 // Plot the histograms
221 //
222 if(ifig >= fNRefFigures) return kFALSE;
770382d9 223 if(!gPad) return kFALSE;
224 gPad->SetLogx(kTRUE);
7fe4e88b 225 TH1 *h(NULL);
1ee39b3a 226 if(ifig < 2){
7fe4e88b 227 if(!(h = (dynamic_cast<TH1 *>(fContainer->At(ifig))))) return kFALSE;
228 h->Draw("e1");
1ee39b3a 229 return kTRUE;
230 }
4226db3e 231 TLegend *leg=new TLegend(.65, .12, .85, .3);
232 leg->SetHeader("Charge");
233 leg->SetBorderSize(1);leg->SetFillColor(kWhite);
1ee39b3a 234 switch(ifig){
235 case 2:
7fe4e88b 236 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram)))) return kFALSE;;
4226db3e 237 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 238 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+1)))) return kFALSE;;
4226db3e 239 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
240 leg->Draw();
1ee39b3a 241 break;
242 case 3:
7fe4e88b 243 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+2)))) return kFALSE;;
4226db3e 244 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 245 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+3)))) return kFALSE;;
4226db3e 246 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
247 leg->Draw();
1ee39b3a 248 break;
249 case 4:
7fe4e88b 250 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+4)))) return kFALSE;;
4226db3e 251 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 252 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+5)))) return kFALSE;;
4226db3e 253 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
254 leg->Draw();
1ee39b3a 255 break;
256 case 5:
7fe4e88b 257 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+6)))) return kFALSE;;
4226db3e 258 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 259 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+7)))) return kFALSE;;
4226db3e 260 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
261 leg->Draw();
770382d9 262 break;
263 case 6:
7fe4e88b 264 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+8)))) return kFALSE;;
4226db3e 265 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 266 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+9)))) return kFALSE;;
4226db3e 267 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
268 leg->Draw();
770382d9 269 break;
270 case 7:
7fe4e88b 271 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+10)))) return kFALSE;;
4226db3e 272 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 273 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+11)))) return kFALSE;;
4226db3e 274 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
275 leg->Draw();
770382d9 276 break;
277 case 8:
7fe4e88b 278 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+12)))) return kFALSE;;
4226db3e 279 h->Draw("e1"); leg->AddEntry(h, " -", "pl");
7fe4e88b 280 if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+13)))) return kFALSE;;
4226db3e 281 h->Draw("e1same"); leg->AddEntry(h, " +", "pl");
282 leg->Draw();
1ee39b3a 283 break;
284 }
285 return kTRUE;
286}
287
288//_____________________________________________________________________________
289TObjArray *AliTRDefficiencyMC::Histos(){
290 //
291 // Create the histograms
292 //
1ee39b3a 293
294 if(fContainer) return fContainer;
770382d9 295 const Int_t nbins = AliTRDCalPID::kNMom;
296 Float_t xbins[nbins+1] = {fgPCut, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
297 const Int_t marker[2][AliPID::kSPECIES+1] = {
298 {20, 21, 22, 23, 29, 2},
299 {24, 25, 26, 27, 30, 5}
300 };
301
302 fContainer = new TObjArray();fContainer->Expand(14);
303
304 TH1 *h(NULL);
305 fContainer->AddAt(h=new TProfile("hEff", "Tracking Efficiency ALL", nbins, xbins), kEfficiencyHistogram);
306 h->SetMarkerStyle(22);
307 h->SetMarkerColor(kBlue);
308 h->GetXaxis()->SetTitle("p [GeV/c]");
309 h->GetXaxis()->SetMoreLogLabels();
310 h->GetYaxis()->SetTitle("Efficiency");
311 h->GetYaxis()->SetRangeUser(0.2, 1.1);
312 fContainer->AddAt(h=new TProfile("hFake", "Fake Tracks", nbins, xbins), kContaminationHistogram);
313 h->SetMarkerStyle(22);
314 h->SetMarkerColor(kBlue);
315 h->GetXaxis()->SetTitle("p [GeV/c]");
316 h->GetXaxis()->SetMoreLogLabels();
317 h->GetYaxis()->SetTitle("Contamination");
318
319 Char_t sign[]={'+', '-'};
320 for(Int_t isign = 0; isign < 2; isign++){
321 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
322 fContainer->AddAt(h=new TProfile(
323 Form("hEff_%s%c", AliPID::ParticleShortName(ispec), sign[isign]),
4226db3e 324 Form("Tracking Efficiency for %s", AliPID::ParticleName(ispec)), nbins, xbins),
770382d9 325 kEfficiencySpeciesHistogram+ispec*2+isign);
326 h->SetMarkerStyle(marker[isign][ispec]);
327 h->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
328 h->SetMarkerColor(kBlack);
329 h->GetXaxis()->SetTitle("p [GeV/c]");
330 h->GetXaxis()->SetMoreLogLabels();
331 h->GetYaxis()->SetTitle("Efficiency");
332 h->GetYaxis()->SetRangeUser(0.2, 1.1);
333 }
334
4226db3e 335 fContainer->AddAt(h=new TProfile(Form("hEff_PID%c", sign[isign]), "Tracking Efficiency no PID", nbins, xbins), kEfficiencySpeciesHistogram+AliPID::kSPECIES*2+isign);
770382d9 336 h->SetMarkerStyle(marker[isign][AliPID::kSPECIES]);
337 h->SetMarkerColor(kBlack);h->SetLineColor(kBlack);
338 h->GetXaxis()->SetTitle("p [GeV/c]");
339 h->GetXaxis()->SetMoreLogLabels();
340 h->GetYaxis()->SetTitle("Efficiency");
341 h->GetYaxis()->SetRangeUser(0.2, 1.1);
342 }
1ee39b3a 343 return fContainer;
344}
345
346//_____________________________________________________________________________
770382d9 347Int_t AliTRDefficiencyMC::IsFindableNot(AliTRDtrackInfo * const trkInf){
1ee39b3a 348 //
349 // Apply Cuts on the Monte Carlo track references
350 // return whether track is findable or not
351 //
1ee39b3a 352
770382d9 353
354 const Float_t chmbHght = AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
355 const Float_t eps(1.E-3);
356 Int_t ntr(trkInf->GetNTrackRefs());
357
358 AliDebug(10, Form(" CANDIDATE TrackRefs[%d]", ntr));
1ee39b3a 359 // Check if track is findable
770382d9 360 Double_t mom(0.), phi(0.), tht(0.);
1ee39b3a 361 Float_t xmin = 10000.0, xmax = 0.0;
362 Float_t ymin = 0.0, ymax = 0.0;
363 Float_t zmin = 0.0, zmax = 0.0;
364 Float_t lastx = 0.0, x = 0.0;
770382d9 365 Int_t nLayers(0), ntrTRD(0);
1ee39b3a 366 Int_t sector[20];
770382d9 367 AliTrackReference *trackRef(NULL);
368 for(Int_t itr(0); itr<ntr; itr++){
369 if(!(trackRef = trkInf->GetTrackRef(itr))) continue;
1ee39b3a 370 x = trackRef->LocalX();
1ee39b3a 371 // Be Sure that we are inside TRD
3d2a3dff 372 if(x < AliTRDinfoGen::GetEndTPC() || x > AliTRDinfoGen::GetEndTRD()) continue;
770382d9 373 sector[ntrTRD] = Int_t(trackRef->Alpha()/AliTRDgeometry::GetAlpha());
374 AliDebug(10, Form(" [%2d] x[%7.2f] y[%7.2f] z[%7.2f] Sec[%2d]", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z(), sector[ntrTRD]));
1ee39b3a 375 if(x < xmin){
376 xmin = trackRef->LocalX();
377 ymin = trackRef->LocalY();
378 zmin = trackRef->Z();
770382d9 379 mom = trackRef->P();
1ee39b3a 380 } else if(x > xmax){
381 xmax = trackRef->LocalX();
382 ymax = trackRef->LocalY();
383 zmax = trackRef->Z();
384 }
385 if(itr > 0){
386 Float_t dist = TMath::Abs(x - lastx);
770382d9 387 if(TMath::Abs(dist - chmbHght) < eps && sector[ntrTRD]==sector[0]){
388 AliDebug(10, Form(" dx = %7.2f", dist));
389 nLayers++;
390 }
1ee39b3a 391 }
392 lastx = x;
770382d9 393 ntrTRD++; if(ntrTRD>=20) break;
394 }
395 Double_t dx(xmax - xmin);
396 if(TMath::Abs(dx)<eps) return kNoChmb;
397
398 phi = (ymax -ymin)/dx;
399 tht = (zmax -zmin)/dx;
400 phi=TMath::ATan(phi)*TMath::RadToDeg();
401 tht=TMath::ATan(tht)*TMath::RadToDeg();
402 Bool_t primary = trkInf->IsPrimary();
403 const AliTRDtrackInfo::AliESDinfo *esd(trkInf->GetESDinfo());
404 AliDebug(10, Form(" p=%6.3f[GeV/c] phi=%6.2f[deg] theta=%6.2f[deg] nLy[%d]",
405 mom, phi, tht, nLayers));
406 if(DebugLevel()){
407 (*DebugStream()) << "IsFindable"
408 << "P=" << mom
409 << "Phi=" << phi
410 << "Tht=" << tht
411 << "Ntr=" << ntrTRD
412 << "NLy=" << nLayers
413 << "Primary=" << primary
414 << "\n";
1ee39b3a 415 }
416
417 // Apply cuts
770382d9 418 if(!nLayers) return kNoChmb;
419 if(xmax < xmin) return kCurved;
420 if(mom < fgPCut) return kPCut;
421
422
423 if(TMath::Abs(phi) > fgPhiCut) return kPhiCut;
424 if(TMath::Abs(tht) > fgThtCut) return kThtCut;
425
426 if(nLayers < 4){
427 if(!esd)return kLayer;
428 if(!(esd->GetStatus() & AliESDtrack::kTPCout)) return kLayer;
1ee39b3a 429 }
770382d9 430
431 //if(!trkInf->IsPrimary()) {failCode=kPrimary; return kFALSE;}
432
433 return kFindable;
1ee39b3a 434}
435
436//_____________________________________________________________________________
770382d9 437void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, ETRDefficiencyMCstatus mode){
1ee39b3a 438 //
439 // Fill Histograms in three different modes:
440 // 1st tracks which are found and accepted
441 // 2nd tracks which are not found and not already counted
442 // 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
443 //
770382d9 444
445 TDatabasePDG *dbPDG(TDatabasePDG::Instance());
446 Double_t trkmom(0.); // the track momentum
447 Int_t trkpdg(-1); // particle PDG code
448 AliTRDtrackInfo *trkInf(NULL);
1ee39b3a 449 for(Int_t itk = 0; itk < nTracks; itk++){
7fe4e88b 450 if(!(trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk])))) continue;
1ee39b3a 451 if(trkInf->GetNTrackRefs()){
452 // use Monte-Carlo Information for Momentum and PID
453 trkmom = trkInf->GetTrackRef(0)->P();
770382d9 454 trkpdg = trkInf->GetPDG();
1ee39b3a 455 }else{
456 // Use TPC Momentum
457 trkmom = trkInf->GetTrack()->P();
458 }
770382d9 459
7fe4e88b 460 TProfile *hp(NULL);
770382d9 461 const Char_t *cmode(NULL);
1ee39b3a 462 switch(mode){
770382d9 463 case kAccept:
7fe4e88b 464 if(!(hp = (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram))))) continue;
465 hp->Fill(trkmom, 1);
466 if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
467 hp->Fill(trkmom, 0);
770382d9 468 cmode="ACCEPT";
1ee39b3a 469 break;
770382d9 470 case kMiss:
7fe4e88b 471 if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))))) continue;
472 hp->Fill(trkmom, 0);
473 if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
474 hp->Fill(trkmom, 0);
770382d9 475 cmode="MISS";
1ee39b3a 476 break;
770382d9 477 case kFake:
7fe4e88b 478 if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
479 hp->Fill(trkmom, 1);
770382d9 480 cmode="FAKE";
1ee39b3a 481 break;
482 }
770382d9 483 AliDebug(3, Form(" track[%d] MC[%d] Mode[%s]", indices[itk], trkInf->GetLabel(), cmode));
484
1ee39b3a 485 // Fill species histogram
770382d9 486 Int_t idxSpec = AliTRDpidUtil::Pdg2Pid(TMath::Abs(trkpdg));
487 Int_t sign = dbPDG->GetParticle(trkpdg)->Charge() > 0. ? 1 : 0;
488 //printf("[%d]%s pdg[%d] sign[%d]\n", idxSpec, AliPID::ParticleName(idxSpec), trkpdg, sign);
489 if(idxSpec < 0) idxSpec = AliPID::kSPECIES;
7fe4e88b 490 if(!(hp = (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + idxSpec*2+sign))))) continue;
491 hp->Fill(trkmom, mode==kAccept?1:0);
1ee39b3a 492 }
493}
494
495//_____________________________________________________________________________
496void AliTRDefficiencyMC::FillStreamTrackWOMC(AliTRDtrackInfo * const trkInf){
497 // fill debug stream
498 // we want to know:
499 // 1. The event number
500 // 2. The track label
501 // 3. The TRD track label
502 // 4. The frequency of the TRD Label
503 // 5. Momentum from TPC (NO STAND ALONE TRACK)
504 // 6. TPC Phi angle
505 // 7. the TRD track
506 // 8. Monte Carlo PID
507 // 9. We check the Labels of the TRD track according to them we search the maching Monte-Carlo track.
508 // From the matching Monte-Carlo track we store trackRefs, phi and momentum
509 // 10. We may also want to keep the kink index
510 Double_t mom = trkInf->GetESDinfo()->GetOuterParam()->P();
511 Int_t event = (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry();
512 Int_t label = trkInf->GetLabel();
513 Int_t kinkIndex = trkInf->GetKinkIndex();
514 Int_t pdg = trkInf->GetPDG();
515 Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
516 Int_t labelsTRD[180]; // Container for the cluster labels
517 Int_t sortlabels[360]; // Cluster Labels sorted according their occurancy
770382d9 518 AliTRDseedV1 *tracklet(NULL);
519 AliTRDcluster *c(NULL);
520 Int_t nclusters(0);
1ee39b3a 521 AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
522 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
523 tracklet = trackTRD->GetTracklet(il);
524 if(!tracklet) continue;
525 tracklet->ResetClusterIter();
770382d9 526 c = NULL;
1ee39b3a 527 while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
528 }
529 // Determine Label and Frequency
530 AliMathBase::Freq(nclusters, const_cast<const Int_t *>(&labelsTRD[0]), &sortlabels[0], kTRUE);
531 Int_t labelTRD = sortlabels[0];
532 Int_t freqTRD = sortlabels[1];
533 // find the track info object matching to the TRD track
534 AliTRDtrackInfo *realtrack = 0;
535 TObjArrayIter rtiter(fTracks);
536 while((realtrack = (AliTRDtrackInfo *)rtiter())){
537 if(realtrack->GetLabel() != labelTRD) continue;
538 break;
539 }
540 TClonesArray trackRefs("AliTrackReference");
541 Int_t realPdg = -1;
542 Double_t realP = 0.;
543 Double_t realPhi = 0.;
544 if(realtrack){
545 // pack the track references into the trackRefsContainer
546 for(Int_t iref = 0; iref < realtrack->GetNTrackRefs(); iref++){
547 new(trackRefs[iref])AliTrackReference(*(realtrack->GetTrackRef(iref)));
548 }
549 realPdg = realtrack->GetPDG();
550 if(realtrack->GetNTrackRefs()){
551 realP = realtrack->GetTrackRef(0)->P();
552 realPhi = realtrack->GetTrackRef(0)->Phi();
553 }
554 }
770382d9 555 (*DebugStream()) << "EffMCfake"
1ee39b3a 556 << "Event=" << event
557 << "Label=" << label
558 << "labelTRD=" << labelTRD
559 << "FreqTRDlabel=" << freqTRD
560 << "TPCp=" << mom
561 << "phiTPC=" << phiTPC
562 << "trackTRD=" << trackTRD
563 << "PDG=" << pdg
564 << "TrackRefs=" << &trackRefs
565 << "RealPDG=" << realPdg
566 << "RealP=" << realP
567 << "RealPhi" << realPhi
568 << "KinkIndex=" << kinkIndex
569 << "\n";
570}
571
572//_____________________________________________________________________________
64d57299 573Int_t AliTRDefficiencyMC::IsRegistered(const AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
1ee39b3a 574 //
575 // Checks if track is registered in a given mode
576 //
770382d9 577
7fe4e88b 578 AliTRDtrackInfo *ti(NULL);
770382d9 579 Int_t label(trkInf->GetLabel());
580 for(Int_t il(nTracks); il--;){
7fe4e88b 581 if(!(ti = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))) continue;
582 if(ti->GetLabel() == label) return il;
1ee39b3a 583 }
770382d9 584 return -1;
1ee39b3a 585}
586