]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/testITSU/AliTaskITSUPerf.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing printout (Davide)
[u/mrichter/AliRoot.git] / ITS / UPGRADE / testITSU / AliTaskITSUPerf.cxx
CommitLineData
fa10255f 1/*************************************************************************
2* Copyright(c) 1998-2008, 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///////////////////////////////////////////////////////////////////////////
17// Class AliTaskITSUPerf //
18// Analysis task to produce data and MC histos needed for tracklets //
19// dNdEta extraction in multiple bins in one go //
20// Author: ruben.shahoyan@cern.ch //
21///////////////////////////////////////////////////////////////////////////
22
23#include "TFile.h"
24#include "TChain.h"
25#include "TTree.h"
26#include "TRandom.h"
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TList.h"
30#include "TObjArray.h"
31#include "TGeoGlobalMagField.h"
32
33#include "AliAnalysisManager.h"
34
35#include "AliESDEvent.h"
36#include "AliESDInputHandler.h"
37#include "AliESDInputHandlerRP.h"
38#include "AliCDBPath.h"
39#include "AliCDBManager.h"
40#include "AliCDBEntry.h"
41#include "AliCDBStorage.h"
42#include "AliGeomManager.h"
43#include "AliMagF.h"
44#include "AliMCEventHandler.h"
45#include "AliMCEvent.h"
46#include "AliMCParticle.h"
47#include "AliStack.h"
48#include "AliGenEventHeader.h"
49#include "AliCentrality.h"
50#include "AliLog.h"
51#include "AliPhysicsSelection.h"
52#include "AliGenEventHeader.h"
53#include "AliGenHijingEventHeader.h"
54#include "AliGenDPMjetEventHeader.h"
55#include "AliGenCocktailEventHeader.h"
56#include "AliESDtrackCuts.h"
57
58#include "AliTaskITSUPerf.h"
59#include "../ITS/UPGRADE/AliITSURecoDet.h"
60#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
61#include "../ITS/UPGRADE/AliITSUClusterPix.h"
62#include "../ITS/UPGRADE/AliITSUTrackCond.h"
7bc38dea 63#include "../ITS/UPGRADE/AliITSUAux.h"
fa10255f 64#include <TGeoManager.h>
7bc38dea 65#include <TTree.h>
fa10255f 66
67
68ClassImp(AliTaskITSUPerf)
69
70
71const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = {
72 "ITSokTPCok"
73 ,"ITSokTPCfk"
74 ,"ITSfkTPCok"
75 ,"ITSfkTPCfk"
76 ,"ITSTPCmismatch"
77 ,"ITSTPCnoMatch"
78};
79
7bc38dea 80
81typedef struct {
82 Bool_t prim;
83 Bool_t rcbl;
84 Char_t nClITS;
85 Char_t nClTPC;
be806700 86 Char_t nClITSMC;
7bc38dea 87 Char_t mcCl[7];
88 Char_t rcCl[7];
89 Char_t fcCl[7];
be806700 90 Char_t charge;
b824a5f5 91 Float_t zV;
7bc38dea 92 Float_t ptMC;
93 Float_t etaMC;
94 Float_t pt;
95 Float_t eta;
b824a5f5 96 Float_t ptTPC;
97 Float_t etaTPC;
7bc38dea 98 Float_t dcaR;
99 Float_t dcaZ;
100 Float_t dcaRE;
101 Float_t dcaZE;
699d8344 102 Float_t ptE;
b824a5f5 103 Float_t ptTPCE;
7bc38dea 104 Float_t phi;
105 Int_t pdg;
106 Int_t lbl;
b824a5f5 107 Int_t lblTPC;
be806700 108 Int_t spl;
b824a5f5 109 Int_t evID;
7bc38dea 110} trInfo_t;
111
112trInfo_t trackInfo;
113
fa10255f 114//________________________________________________________________________
115/*//Default constructor
116AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
117 : AliAnalysisTaskSE(name),
118*/
119//________________________________________________________________________
120AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
121 : AliAnalysisTaskSE(name)
122 //
123 ,fOutput(0)
124 ,fHistosCentMCLb(0)
125 ,fHistosCent(0)
126 ,fRPTree(0)
127 ,fStack(0)
128 ,fMCEvent(0)
129 ,fVtxSPD(0)
130 ,fVtxTrc(0)
131 ,fESDEvent(0)
132 ,fUseSpecialOutput(kTRUE)
133 ,fUseMC(kTRUE)
134 ,fTrackingCond(0)
135 //
136 ,fGeom(0)
137 ,fITS(0)
138 ,fNSelTracksMC(0)
139 ,fMCStatus(0)
097e4cdf 140 //
141#ifdef _CLUS_LIST_
142 ,fClLists(0)
143#endif
144 //
fa10255f 145 ,fNPtBins(20)
146 ,fNResBins(100)
097e4cdf 147 ,fMinTPCclusters(60)
fa10255f 148 ,fPtMin(0)
149 ,fPtMax(10.)
150 ,fEtaMin(-3.)
151 ,fEtaMax( 3.)
152 ,fZVertexMin(-15.)
153 ,fZVertexMax(15.)
154 //
155 ,fCurrCentBin(0)
156 ,fNCentBins(1)
157 ,fUseCentralityVar(0)
60b14b59 158 ,fTPCCut(0)
7bc38dea 159 ,fTree(0)
fa10255f 160{
161 // Constructor
162
163 DefineOutput(1, TList::Class());
164 //
165}
166
167//________________________________________________________________________
168AliTaskITSUPerf::~AliTaskITSUPerf()
169{
170 // Destructor
171 // histograms are in the output list and deleted when the output
172 // list is deleted by the TSelector dtor
173 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
174 printf("Deleteing output\n");
175 delete fOutput;
176 fOutput = 0;
177 }
178 //
179 fHistosCentMCLb.Clear();
180 fHistosCent.Clear();
181 //
182 delete fITS;
183 delete fGeom;
60b14b59 184 delete fTPCCut;
fa10255f 185 //
186}
187
188//________________________________________________________________________
189void AliTaskITSUPerf::UserCreateOutputObjects()
190{
191 //
192 if (fUseSpecialOutput) OpenFile(1);
193 fOutput = new TList();
194 fOutput->SetOwner();
195 //
196 // init geometry
197 if (!gGeoManager) {
198 TString geom = "geometry.root";
199 AliGeomManager::LoadGeometry(geom.Data());
200 if (!gGeoManager) AliFatal("Failed to load geometry");
201 }
202 //
203 // create ITS interface
204 fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE);
205 AliITSUClusterPix::SetGeom(fGeom);
206 fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface");
207 fITS->CreateClusterArrays();
208 //
209 // Create histograms
210 for (int ib=0;ib<fNCentBins;ib++) {
211 BookHistos(ib);
212 }
213 //
214 int nhist = fOutput->GetEntries();
215 for (int i=0;i<nhist;i++) {
216 TObject* hst = fOutput->At(i);
217 if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
218 ((TH1*)hst)->Sumw2();
219 }
220 //
7bc38dea 221 fTree = new TTree("trinfo","trinfo");
222 fTree->Branch("prim",&trackInfo.prim,"prim/O");
223 fTree->Branch("rcbl",&trackInfo.rcbl,"rcbl/O");
224 fTree->Branch("nClITS",&trackInfo.nClITS,"nClITS/b");
225 fTree->Branch("nClTPC",&trackInfo.nClTPC,"nClTPC/b");
be806700 226 fTree->Branch("nClITSMC",&trackInfo.nClITSMC,"nClITSMC/b");
7bc38dea 227 fTree->Branch("mcCl",&trackInfo.mcCl,"mcCl[7]/b");
228 fTree->Branch("rcCl",&trackInfo.rcCl,"rcCl[7]/b");
229 fTree->Branch("fcCl",&trackInfo.fcCl,"fcCl[7]/b");
be806700 230 fTree->Branch("charge",&trackInfo.charge,"charge/B");
b824a5f5 231 fTree->Branch("zV", &trackInfo.zV,"zV/F");
7bc38dea 232 fTree->Branch("ptMC", &trackInfo.ptMC,"ptMC/F");
233 fTree->Branch("etaMC", &trackInfo.etaMC,"etaMC/F");
234 fTree->Branch("pt", &trackInfo.pt,"pt/F");
235 fTree->Branch("eta", &trackInfo.eta,"eta/F");
b824a5f5 236 fTree->Branch("ptTPC", &trackInfo.ptTPC,"ptTPC/F");
237 fTree->Branch("etaTPC", &trackInfo.etaTPC,"etaTPC/F");
7bc38dea 238 fTree->Branch("dcaR", &trackInfo.dcaR,"dcaR/F");
239 fTree->Branch("dcaZ", &trackInfo.dcaZ,"dcaZ/F");
240 fTree->Branch("dcaRE", &trackInfo.dcaRE,"dcaRE/F");
241 fTree->Branch("dcaZE", &trackInfo.dcaZE,"dcaZE/F");
699d8344 242 fTree->Branch("ptE", &trackInfo.ptE,"ptE/F");
b824a5f5 243 fTree->Branch("ptTPCE", &trackInfo.ptTPCE,"ptTPCE/F");
7bc38dea 244 fTree->Branch("phi", &trackInfo.phi,"phi/F");
b824a5f5 245 fTree->Branch("pdg", &trackInfo.pdg,"pdg/I");
7bc38dea 246 fTree->Branch("lbl", &trackInfo.lbl,"lbl/I");
b824a5f5 247 fTree->Branch("lblTPC", &trackInfo.lblTPC,"lblTPC/I");
be806700 248 fTree->Branch("spl", &trackInfo.spl,"spl/I");
b824a5f5 249 fTree->Branch("evID", &trackInfo.evID,"evID/I");
7bc38dea 250 fOutput->Add(fTree);
251 //
60b14b59 252 fTPCCut = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
253 fTPCCut->SetEtaRange(fEtaMin,fEtaMax);
254 fTPCCut->SetPtRange(fPtMin,fPtMax);
255 fTPCCut->SetPtRange(fPtMin,fPtMax);
256 fTPCCut->SetMinNClustersTPC(fMinTPCclusters);
7bc38dea 257 printf("Cuts: Eta: %f %f | Pt: %f %f\n",fEtaMin,fEtaMax,fPtMin,fPtMax);
60b14b59 258 //
097e4cdf 259#ifdef _CLUS_LIST_
260 fClLists.SetOwner();
261#endif
fa10255f 262 PostData(1, fOutput);
263 //
264}
265
266//________________________________________________________________________
267void AliTaskITSUPerf::UserExec(Option_t *)
268{
269 // Main loop
270 static int evID = 0;
271 AliInfo(Form("Event: %d",evID++));
272 //
273 AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
274 fRPTree = 0;
275 AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
b824a5f5 276 if (!handRP) { AliFatal("No RP handler"); return; }
fa10255f 277 //
278 fESDEvent = handRP->GetEvent();
279 if (!fESDEvent) { AliFatal("No AliESDEvent"); return; }
280 //
281 fRPTree = handRP->GetTreeR("ITS");
b824a5f5 282 if (!fRPTree) { AliWarning("Invalid ITS cluster tree"); }
fa10255f 283 //
284 AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
285 if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; }
286 fMCEvent = eventHandler->MCEvent();
287 if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; }
288 fStack = fMCEvent->Stack();
289 if (!fStack) { AliFatal("Stack is not available"); return; }
290 //
291 // MC Generator info
292 AliGenEventHeader* mcGenH = 0;
293 mcGenH = fMCEvent->GenEventHeader();
294 TArrayF vtmc(3);
295 mcGenH->PrimaryVertex(vtmc);
296 for (int i=3;i--;) fVtxMC[i] = vtmc[i];
297 //
298 printf("MagField: %f\n",fESDEvent->GetMagneticField());
299 printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]);
300 //
301 fVtxSPD = fESDEvent->GetPrimaryVertexSPD();
302 //
303 if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e | Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors());
304 else printf("VtxSPD: N/A\n");
305 //
306 fVtxTrc = fESDEvent->GetPrimaryVertexTracks();
307 if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e | Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors());
308 else printf("VtxTRC: N/A\n");
309 //
310 BuildMCInfo();
311 CheckTracks();
312 //
313}
314
315
316//________________________________________________________________________
317void AliTaskITSUPerf::Terminate(Option_t *)
318{
097e4cdf 319#ifdef _CLUS_LIST_
320 fClLists.Delete();
321#endif
fa10255f 322 Printf("Terminating...");
323}
324
325//_________________________________________________________________________
326void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb)
327{
328 // book standard set of histos for specific centrality bin and MC status type,
329 // adding them to output list and management array
330 //
60b14b59 331 const double kMaxDPt = 0.2;
332 const double kMaxDCA = 0.1;
fa10255f 333
334 TH2* h2=0;
335 // pt resolution
336 h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"),
337 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"),
338 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt);
339 AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) );
340 //
341 // transverse DCA resolution
342 h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"),
343 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"),
344 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
345 AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) );
346 //
347 // Z DCA resolution
348 h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"),
349 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"),
350 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
351 AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) );
352 //
353}
354
355//_________________________________________________________________________
356void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin)
357{
358 // book standard set of histos for specific centrality bin
359 //
360 TH2* h2=0;
361 //
362 // MC labels combination vs pt
097e4cdf 363 h2 = new TH2F(Form("B%d_%s_RCBL",bin,"MCLabvsPT"),
364 Form("B%d %s Reconstructable",bin,"MCLabComb vs p_{T}"),
fa10255f 365 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
097e4cdf 366 AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusRcbl,-1,bin) );
fa10255f 367 for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
368 //
097e4cdf 369 h2 = new TH2F(Form("B%d_%s_NotRCBL_prim",bin,"MCLabvsPT"),
370 Form("B%d %s Not Reconstructable, Prim.",bin,"MCLabComb vs p_{T}"),
371 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
372 AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblPrim,-1,bin) );
373 for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
374 //
375 h2 = new TH2F(Form("B%d_%s_NotRCBL_sec",bin,"MCLabvsPT"),
376 Form("B%d %s Not Reconstructable, Sec.",bin,"MCLabComb vs p_{T}"),
377 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
378 AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblSec,-1,bin) );
379 for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
7bc38dea 380 //
381 h2 = new TH2F(Form("B%d_MCLr_Prim",bin),
382 Form("B%d MCLr Prim vs pt",bin),
383 fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5);
384 AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresPrim,-1,bin) );
385 //
386 h2 = new TH2F(Form("B%d_MCLr_Sec",bin),
387 Form("B%d MCLr Sec vs pt",bin),
388 fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5);
389 AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresSec,-1,bin) );
390 //
fa10255f 391}
392
393//_________________________________________________________________________
394void AliTaskITSUPerf::BookHistos(Int_t bin)
395{
396 // book standard set of histos, adding them to output list and management array
397 //
398 // standard histos
399 for (int ilb=0;ilb<kNLabelTypes;ilb++) { // create similar set of histos for all kind of MC labels
400 BookStandardHistosCentMCLb(bin,ilb);
401 }
402 BookStandardHistosCent(bin); // st.histos regardless of MClabels comb.
403}
404
405//_________________________________________________________________________
406void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at)
407{
408 // add single histo to the set
409 if (at>=0) array->AddAtAndExpand(h,at);
410 else array->Add(h);
411 fOutput->Add(h);
412}
413
414//_________________________________________________________________________
415Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const
416{
417 // generate ID for the histo. htype must be <100, mcStat<kNLabelTypes
418 //
7bc38dea 419 // if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes));
fa10255f 420 if (centBin>=fNCentBins) AliFatal(Form("CentBin ID=%d > max.allowed %d",centBin,fNCentBins));
421 //
422 if (centBin>=0) {
423 if (mcStat>=0) { // MCtruthness dependent histos
424 if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC));
425 return htype + kHNStdHistosCentMC*(mcStat + kNLabelTypes*centBin);
426 }
427 else {
428 if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent));
429 return htype + kHNStdHistosCent*centBin;
430 }
431 }
432 // custom histo
433 return htype;
434 //
435}
436
437//_________________________________________________________________________
438void AliTaskITSUPerf::CheckTracks()
439{
440 // build mc truth info for tracks
097e4cdf 441 Bool_t dump = kFALSE;
7bc38dea 442 Bool_t dump1 = kFALSE;
b824a5f5 443 static int evid = 0;
444 trackInfo.evID = evid++;
445 trackInfo.zV = fESDEvent->GetPrimaryVertex()->GetZ();
fa10255f 446 //
447 int ntr = fESDEvent->GetNumberOfTracks();
448 int ntrMC = fStack->GetNtrack();
449 //
450 double fldBz = fESDEvent->GetMagneticField();
451 float dcaRZ[2];
452 //
453 for (int itr=0;itr<ntr;itr++) {
454 AliESDtrack* trc = fESDEvent->GetTrack(itr);
455 //
456 // at the moment we consider only TPC/ITS tracks
60b14b59 457 //if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
458 if (!fTPCCut->IsSelected(trc)) continue;
7bc38dea 459 //
460 trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ);
461 //
fa10255f 462 Int_t labMC = trc->GetLabel();
463 Int_t labMCabs = TMath::Abs(labMC);
464 Int_t labMCTPC = trc->GetTPCLabel();
465 Int_t labMCITS = trc->GetITSLabel();
466 Int_t nClTPC = trc->GetNcls(1);
467 Int_t nClITS = trc->GetNcls(0);
468 //
469 UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs];
097e4cdf 470 double pt = trc->Pt();
471 double eta = trc->Eta();
472 //
7bc38dea 473 TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec));
474 if (hlr) {
475 for (int il=0;il<7;il++) {
476 if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(pt,il);
477 }
478 hlr->Fill(pt,-99); // count tracks
479 }
480 //
481
fa10255f 482 Bool_t reject = mcStatus & BIT(kTrCondFail);
fa10255f 483 //
484 int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS);
485 //
60b14b59 486 // if (eta>fEtaMax || eta<fEtaMin) continue;
487 // if (nClTPC<fMinTPCclusters) continue;
7bc38dea 488 // printf("#%3d pt:%5.2f eta:%+5.2f | Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d | Ntpc:%3d Nits:%3d | Rej:%d\n",
489 // itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS,reject);
490 AliMCParticle *part = 0;
491 if (labMCabs>=ntrMC) continue;
492 part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs);
493 double ptMC = part->Pt();
494 double etaMC = part->Eta();
fa10255f 495 //
7bc38dea 496 if (fTree) {
497 trackInfo.lbl = labMC;
b824a5f5 498 trackInfo.lblTPC = labMCTPC;
7bc38dea 499 trackInfo.prim = (mcStatus&BIT(kMCPrimBit)) != 0;
500 trackInfo.rcbl = !reject;
501 trackInfo.nClITS = trc->GetNcls(0);
502 trackInfo.nClTPC = trc->GetNcls(1);
503 trackInfo.ptMC = ptMC;
504 trackInfo.etaMC = etaMC;
be806700 505 trackInfo.charge = trc->Charge();
7bc38dea 506 trackInfo.pt = trc->Pt();
be806700 507 trackInfo.eta = trc->Eta();
b824a5f5 508 //
509
510 //
be806700 511 trackInfo.spl = trc->GetITSModuleIndex(10);
7bc38dea 512 trackInfo.phi = part->Phi();
513 trackInfo.pdg = part->PdgCode();
514 trackInfo.dcaR = dcaRZ[0];
515 trackInfo.dcaZ = dcaRZ[1];
516 trackInfo.dcaRE = TMath::Sqrt(trc->GetSigmaY2());
517 trackInfo.dcaZE = TMath::Sqrt(trc->GetSigmaZ2());
b824a5f5 518 trackInfo.ptE = TMath::Sqrt(trc->GetSigma1Pt2());
519
520 const AliExternalTrackParam* ttpc = trc->GetTPCInnerParam();
521 if (ttpc) {
522 trackInfo.ptTPC = ttpc->Pt();
523 trackInfo.etaTPC = ttpc->Eta();
524 trackInfo.ptTPCE = TMath::Sqrt(ttpc->GetSigma1Pt2());//*ttpc->Pt()*ttpc->Pt();
525 }
526 else {
527 trackInfo.ptTPC = -1;
528 trackInfo.etaTPC = 0;
529 trackInfo.ptTPCE = -1;
530 }
7bc38dea 531 //
be806700 532 trackInfo.nClITSMC = 0;
7bc38dea 533 for (int il=0;il<7;il++) {
be806700 534 trackInfo.mcCl[il] = (mcStatus & (0x1<<(il+kITSHitBits))) != 0;
535 if (trackInfo.mcCl[il]) trackInfo.nClITSMC++;
b824a5f5 536 trackInfo.rcCl[il] = trc->HasPointOnITSLayer(il);
7bc38dea 537 trackInfo.fcCl[il] = trc->HasSharedPointOnITSLayer(il);
538 }
539 fTree->Fill();
540 }
097e4cdf 541 if (reject) {
542 GetHisto(&fHistosCent,(mcStatus & BIT(kMCPrimBit)) ? kHMatchStatusNRcblPrim:kHMatchStatusNRcblSec)->Fill(pt,mcLabType); // matching status
543 if (dump) {
544 printf("--#%4d Pt:%5.2f Eta:%5.2f |",itr,pt,eta);
545 for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0);
546 printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
547 }
548 continue;
549 }
550 GetHisto(&fHistosCent,kHMatchStatusRcbl)->Fill(pt,mcLabType); // matching status
fa10255f 551 //
097e4cdf 552 TObjArray* clList = 0;
fa10255f 553 //
7bc38dea 554 if (part) {
555 /*
556 int pdg = part->PdgCode();
557 if (TMath::Abs(pdg)!=211) continue;
558 if (nClITS<7) continue;
559 */
097e4cdf 560#ifdef _CLUS_LIST_
097e4cdf 561 //
562 clList = (TObjArray*)fClLists.At(labMCabs);
563 if (clList) {
564 printf("Clusters for track %d MCLabel:%d\n",itr,labMC);
565 for (int ic=0;ic<clList->GetEntriesFast();ic++) {
566 AliITSUClusterPix* cl = (AliITSUClusterPix*) clList->At(ic);
567 cl->Print();
568 }
569 }
097e4cdf 570#endif
7bc38dea 571 }
097e4cdf 572
fa10255f 573 if (part) {
574 //
fa10255f 575 //
097e4cdf 576 if (dump) {
577 printf(" #%4d Pt:%5.2f Eta:%5.2f |",itr,ptMC,etaMC);
60b14b59 578 for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0);
579 printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
097e4cdf 580 }
581 //
582 if ( (mcStatus&BIT(kMCPrimBit)) ) {
be806700 583 if (dump1 && TMath::Abs(ptMC-0.5)<0.1) {
584 printf("#%4d Pt:%5.2f Eta:%5.2f | ",itr,ptMC,etaMC);
585 for (int k=0;k<7;k++) printf("(%d/%d)", (mcStatus&(0x1<<k)) ? 1:0, trc->HasPointOnITSLayer(k));
7bc38dea 586 printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
587 }
fa10255f 588 // compare MC vs reco track params
589 if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC);
590 //
591 // DCA resolution
fa10255f 592 GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]);
593 GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]);
594 //
595 }
596 //
597 }
598 //
599
600
601 }
602 //
603
604}
605
606//_________________________________________________________________________
607void AliTaskITSUPerf::BuildMCInfo()
608{
609 // build mc truth info for tracks
610 //
611 int ntrMC = fStack->GetNtrack();
612 //
613 if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC);
614 fMCStatus.Reset();
097e4cdf 615#ifdef _CLUS_LIST_
616 fClLists.Delete();
617 if (fClLists.GetSize()<ntrMC) fClLists.Expand(ntrMC);
618 //
619 // create empty lists for reconstructed tracks
620 int ntr = fESDEvent->GetNumberOfTracks();
621 for (int itr=0;itr<ntr;itr++) {
622 AliESDtrack* trc = fESDEvent->GetTrack(itr);
623 int lb = TMath::Abs(trc->GetLabel());
624 if (lb<ntrMC && !fClLists.At(itr)) {
625 fClLists.AddAt(new TObjArray(7),lb);
626 }
627 }
628#endif
fa10255f 629 //
630 for (int itr=0;itr<ntrMC;itr++) {
631 //
632 AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(itr);
633 if (!part->Charge()) continue;
634 //
635 UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
636 mcStatus = 0;
637 //
638 if (fStack->IsPhysicalPrimary(itr)) mcStatus |= BIT(kMCPrimBit);
639 //
640 /*
641 TParticle* tp = part->Particle();
642 printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print();
643 */
644 //
645 }
646 //
647 int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0;
648 // fill MC clusters info (if available)
649 if (fRPTree) {
650 // load clusters
651 fITS->LoadClusters(fRPTree);
097e4cdf 652 fITS->ProcessClusters(); // order in the same way as in reconstruction
fa10255f 653 //
097e4cdf 654 for (int ilr=0;ilr<fITS->GetNLayersActive();ilr++) {
fa10255f 655 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
656 for (int icl=lr->GetNClusters();icl--;) {
657 AliCluster* cl = lr->GetCluster(icl);
658 // check labels
659 for (int ilb=0;ilb<3;ilb++) {
660 Int_t lab = cl->GetLabel(ilb);
661 if (lab<0||lab>=ntrMC) continue;
662 UInt_t& mcStatus = (UInt_t &)fMCStatus[lab];
663 mcStatus |= 0x1<<(ilr+kITSHitBits);
097e4cdf 664#ifdef _CLUS_LIST_
665 TObjArray *clList = (TObjArray*)fClLists.At(lab);
666 if (clList) clList->AddLast(cl); // create cl.list only for reconstructed tracks
667#endif
fa10255f 668 }
669 }
670 }
671 //
672 if (nTrCond) {
673 for (int itr=0;itr<ntrMC;itr++) {
674 //
675 UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
676 Bool_t fail = kTRUE;
677 UShort_t patt = (mcStatus&0xffff);
678 for (int ic=0;ic<nTrCond;ic++) {
679 AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic);
680 if (cnd->CheckPattern(patt)) {fail=kFALSE; break;}
681 }
682 if (fail) mcStatus |= BIT(kTrCondFail);
7bc38dea 683 //
684 /*
685 int nlr = AliITSUAux::NumberOfBitsSet(patt);
686 AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(itr);
687 double etaMC = part->Eta();
688 if (etaMC>=fEtaMin && etaMC<=fEtaMax && nlr) {
689 TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec));
690 if (hlr) {
691 double ptMC = part->Pt();
692 for (int il=0;il<7;il++) {
693 if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(ptMC,il);
694 }
695 hlr->Fill(ptMC,-99); // count tracks
696 }
697 }
698 */
fa10255f 699 }
700 }
701 //
702 }
703 //
704}
705
706
707//_________________________________________________________________________
708Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS)
709{
710 // build type of track (0: both ITS and TPC correct, 1: ITS corr. TPC fake, 2: ITS fake TPC corr., 3: ITS fake TPC fake, 4: mismatch
711 //
712 int tp = 0;
713 if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch;
714 if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch;
715 else {
716 if (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1;
717 else tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1;
718 }
719 return tp;
720 //
721}
722
723//_________________________________________________________________________
724Int_t AliTaskITSUPerf::GetCentralityBin() const
725{
726 // calculate centrality bin
727 return 0;
728}
729