]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDinfoGen.cxx
full resolution monitoring for clusters, tracklets and trackIN
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDinfoGen.cxx
CommitLineData
94b94be0 1/**************************************************************************\r
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3* *\r
4* Author: The ALICE Off-line Project. *\r
5* Contributors are mentioned in the code where appropriate. *\r
6* *\r
7* Permission to use, copy, modify and distribute this software and its *\r
8* documentation strictly for non-commercial purposes is hereby granted *\r
9* without fee, provided that the above copyright notice appears in all *\r
10* copies and that both the copyright notice and this permission notice *\r
11* appear in the supporting documentation. The authors make no claims *\r
12* about the suitability of this software for any purpose. It is *\r
13* provided "as is" without express or implied warranty. *\r
14**************************************************************************/\r
15\r
16/* $Id: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */\r
17\r
18////////////////////////////////////////////////////////////////////////////\r
19//\r
20// Tender wagon for TRD performance/calibration train\r
21//\r
22// In this wagon the information from\r
23// - ESD\r
24// - Friends [if available]\r
25// - MC [if available]\r
26// are grouped into AliTRDtrackInfo objects and fed to worker tasks\r
27//\r
28// Authors:\r
29// Markus Fasel <M.Fasel@gsi.de>\r
30// Alexandru Bercuci <A.Bercuci@gsi.de>\r
31//\r
32////////////////////////////////////////////////////////////////////////////\r
33\r
34#include <TClonesArray.h>\r
35#include <TObjArray.h>\r
36#include <TObject.h>\r
37#include <TString.h>\r
38#include <TH1S.h>\r
39#include <TPad.h>\r
40#include <TFile.h>\r
41#include <TTree.h>\r
42#include <TROOT.h>\r
43#include <TChain.h>\r
44#include <TParticle.h>\r
45\r
46#include "AliLog.h"\r
47#include "AliAnalysisManager.h"\r
48#include "AliGeomManager.h"\r
49#include "AliCDBManager.h"\r
4e38ea4c 50#include "AliCDBStorage.h"\r
94b94be0 51#include "AliCDBEntry.h"\r
52#include "AliCDBPath.h"\r
53#include "AliESDEvent.h"\r
54#include "AliMCEvent.h"\r
55#include "AliESDInputHandler.h"\r
56#include "AliMCEventHandler.h"\r
57\r
58#include "AliESDfriend.h"\r
59#include "AliESDfriendTrack.h"\r
60#include "AliESDHeader.h"\r
61#include "AliESDRun.h"\r
62#include "AliESDtrack.h"\r
63#include "AliESDv0.h"\r
64#include "AliESDtrackCuts.h"\r
65#include "AliMCParticle.h"\r
66#include "AliMultiplicity.h"\r
67#include "AliPID.h"\r
68#include "AliStack.h"\r
69#include "AliTrackReference.h"\r
70#include "TTreeStream.h"\r
71\r
72#include <cstdio>\r
73#include <climits>\r
74#include <cstring>\r
75#include <iostream>\r
76#include <memory>\r
77\r
78#include "AliTRDReconstructor.h"\r
79#include "AliTRDrecoParam.h"\r
80#include "AliTRDcalibDB.h"\r
81#include "AliTRDtrackerV1.h"\r
82#include "AliTRDgeometry.h"\r
83#include "AliTRDtrackV1.h"\r
84#include "AliTRDseedV1.h"\r
85#include "AliTRDcluster.h"\r
86#include "AliTRDinfoGen.h"\r
87#include "AliTRDpwg1Helper.h"\r
88#include "info/AliTRDtrackInfo.h"\r
89#include "info/AliTRDeventInfo.h"\r
90#include "info/AliTRDv0Info.h"\r
91#include "info/AliTRDeventCuts.h"\r
92\r
93ClassImp(AliTRDinfoGen)\r
94\r
95const Float_t AliTRDinfoGen::fgkITS = 100.; // to be checked\r
96const Float_t AliTRDinfoGen::fgkTPC = 290.;\r
97const Float_t AliTRDinfoGen::fgkTRD = 365.;\r
98\r
99const Float_t AliTRDinfoGen::fgkEvVertexZ = 15.;\r
100const Int_t AliTRDinfoGen::fgkEvVertexN = 1;\r
101\r
102const Float_t AliTRDinfoGen::fgkTrkDCAxy = 3.;\r
103const Float_t AliTRDinfoGen::fgkTrkDCAz = 10.;\r
104const Int_t AliTRDinfoGen::fgkNclTPC = 70;\r
105const Float_t AliTRDinfoGen::fgkPt = 0.2;\r
106const Float_t AliTRDinfoGen::fgkEta = 0.9;\r
107AliTRDReconstructor* AliTRDinfoGen::fgReconstructor(NULL);\r
108AliTRDgeometry* AliTRDinfoGen::fgGeo(NULL);\r
109//____________________________________________________________________\r
110AliTRDinfoGen::AliTRDinfoGen()\r
111 :AliAnalysisTaskSE()\r
112 ,fEvTrigger(NULL)\r
113 ,fESDev(NULL)\r
114 ,fMCev(NULL)\r
115 ,fEventCut(NULL)\r
116 ,fTrackCut(NULL)\r
117 ,fV0Cut(NULL)\r
118 ,fOCDB("local://$ALICE_ROOT/OCDB")\r
119 ,fTrackInfo(NULL)\r
120 ,fEventInfo(NULL)\r
121 ,fV0Info(NULL)\r
122 ,fTracksBarrel(NULL)\r
123 ,fTracksSA(NULL)\r
124 ,fTracksKink(NULL)\r
125 ,fV0List(NULL)\r
126 ,fContainer(NULL)\r
246fe7f7 127 ,fRecos(NULL)\r
94b94be0 128 ,fDebugStream(NULL)\r
129{\r
130 //\r
131 // Default constructor\r
132 //\r
133 SetNameTitle("TRDinfoGen", "MC-REC TRD-track list generator");\r
134}\r
135\r
136//____________________________________________________________________\r
137AliTRDinfoGen::AliTRDinfoGen(char* name)\r
138 :AliAnalysisTaskSE(name)\r
139 ,fEvTrigger(NULL)\r
140 ,fESDev(NULL)\r
141 ,fMCev(NULL)\r
142 ,fEventCut(NULL)\r
143 ,fTrackCut(NULL)\r
144 ,fV0Cut(NULL)\r
145 ,fOCDB("local://$ALICE_ROOT/OCDB")\r
146 ,fTrackInfo(NULL)\r
147 ,fEventInfo(NULL)\r
148 ,fV0Info(NULL)\r
149 ,fTracksBarrel(NULL)\r
150 ,fTracksSA(NULL)\r
151 ,fTracksKink(NULL)\r
152 ,fV0List(NULL)\r
153 ,fContainer(NULL)\r
246fe7f7 154 ,fRecos(NULL)\r
94b94be0 155 ,fDebugStream(NULL)\r
156{\r
157 //\r
158 // Default constructor\r
159 //\r
160 SetTitle("MC-REC TRD-track list generator");\r
161 DefineOutput(AliTRDpwg1Helper::kTracksBarrel, TObjArray::Class());\r
162 DefineOutput(AliTRDpwg1Helper::kTracksSA, TObjArray::Class());\r
163 DefineOutput(AliTRDpwg1Helper::kTracksKink, TObjArray::Class());\r
164 DefineOutput(AliTRDpwg1Helper::kEventInfo, AliTRDeventInfo::Class());\r
165 DefineOutput(AliTRDpwg1Helper::kV0List, TObjArray::Class());\r
166 DefineOutput(AliTRDpwg1Helper::kMonitor, TObjArray::Class()); // histogram list\r
167}\r
168\r
169//____________________________________________________________________\r
170AliTRDinfoGen::~AliTRDinfoGen()\r
171{\r
172// Destructor\r
173 if(fgGeo) delete fgGeo;\r
174 if(fgReconstructor) delete fgReconstructor;\r
175 if(fDebugStream) delete fDebugStream;\r
176 if(fEvTrigger) delete fEvTrigger;\r
177 if(fV0Cut) delete fV0Cut;\r
178 if(fTrackCut) delete fTrackCut;\r
179 if(fEventCut) delete fEventCut;\r
180 if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;\r
181 if(fEventInfo) delete fEventInfo; fEventInfo = NULL;\r
182 if(fV0Info) delete fV0Info; fV0Info = NULL;\r
183 if(fTracksBarrel){ \r
184 fTracksBarrel->Delete(); delete fTracksBarrel;\r
185 fTracksBarrel = NULL;\r
186 }\r
187 if(fTracksSA){ \r
188 fTracksSA->Delete(); delete fTracksSA;\r
189 fTracksSA = NULL;\r
190 }\r
191 if(fTracksKink){ \r
192 fTracksKink->Delete(); delete fTracksKink;\r
193 fTracksKink = NULL;\r
194 }\r
195 if(fV0List){ \r
196 fV0List->Delete(); \r
197 delete fV0List;\r
198 fV0List = NULL;\r
199 }\r
200 if(fContainer){ \r
201 fContainer->Delete(); \r
202 delete fContainer;\r
203 fContainer = NULL;\r
204 }\r
205}\r
206\r
207//____________________________________________________________________\r
208Bool_t AliTRDinfoGen::GetRefFigure(Int_t)\r
209{\r
210// General graphs for PWG1/TRD train\r
211 if(!gPad){\r
212 AliWarning("Please provide a canvas to draw results.");\r
213 return kFALSE;\r
214 }\r
3ceb45ae 215 fContainer->At(kStatTrk)->Draw("bar");\r
94b94be0 216 return kTRUE;\r
217}\r
218\r
219//____________________________________________________________________\r
220void AliTRDinfoGen::UserCreateOutputObjects()\r
221{ \r
222 //\r
223 // Create Output Containers (TObjectArray containing 1D histograms)\r
224 //\r
225 \r
226 fTrackInfo = new AliTRDtrackInfo();\r
227 fEventInfo = new AliTRDeventInfo();\r
228 fV0Info = new AliTRDv0Info();\r
229 fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);\r
230 fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);\r
231 fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);\r
232 fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);\r
233\r
234 // define general monitor\r
3ceb45ae 235 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);\r
94b94be0 236 TH1 *h=new TH1I("hStat", "Run statistics;Observable;Entries", Int_t(kNObjects), -0.5, Float_t(kNObjects)-0.5);\r
237 TAxis *ax(h->GetXaxis());\r
238 ax->SetBinLabel(Int_t(kTracksESD) + 1, "ESD");\r
239 ax->SetBinLabel(Int_t(kTracksMC) + 1, "MC");\r
240 ax->SetBinLabel(Int_t(kV0) + 1, "V0");\r
241 ax->SetBinLabel(Int_t(kTPC) + 1, "TPC");\r
242 ax->SetBinLabel(Int_t(kTRDin) + 1, "TRDin");\r
243 ax->SetBinLabel(Int_t(kTRDout) + 1, "TRDout");\r
244 ax->SetBinLabel(Int_t(kBarrel) + 1, "Barrel");\r
245 ax->SetBinLabel(Int_t(kBarrelMC) + 1, "BarrelMC");\r
246 ax->SetBinLabel(Int_t(kSA) + 1, "SA");\r
247 ax->SetBinLabel(Int_t(kSAMC) + 1, "SAMC");\r
248 ax->SetBinLabel(Int_t(kKink) + 1, "Kink");\r
249 ax->SetBinLabel(Int_t(kKinkMC) + 1, "KinkMC");\r
250 ax->SetBinLabel(Int_t(kBarrelFriend) + 1, "BFriend");\r
251 ax->SetBinLabel(Int_t(kSAFriend) + 1, "SFriend");\r
3ceb45ae 252 fContainer->AddAt(h, kStatTrk);\r
246fe7f7 253 h=new TH1I("hEv", "Run statistics;Event Class;Entries", 4, -0.5, 3.5);\r
254 ax = h->GetXaxis();\r
255 ax->SetBinLabel(1, "Low");\r
256 ax->SetBinLabel(2, "High");\r
257 ax->SetBinLabel(3, "Cosmic");\r
258 ax->SetBinLabel(4, "Calib");\r
3ceb45ae 259 fContainer->AddAt(h, kEvType);\r
3ed01fbe 260 h=new TH1I("hBC", "TOF Bunch Cross statistics;BC index;Entries", 31, -10.5, 20.5);\r
3ceb45ae 261 fContainer->AddAt(h, kBunchCross);\r
94b94be0 262 PostData(AliTRDpwg1Helper::kMonitor, fContainer);\r
263}\r
264\r
265//____________________________________________________________________\r
266Bool_t AliTRDinfoGen::Load(const Char_t *file, const Char_t *dir, const Char_t *name)\r
267{\r
268// Load data from performance file\r
269\r
270 if(!TFile::Open(file)){\r
271 AliWarning(Form("Couldn't open file %s.", file));\r
272 return kFALSE;\r
273 }\r
274 if(dir){\r
275 if(!gFile->cd(dir)){\r
276 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
277 return kFALSE;\r
278 }\r
279 }\r
280 TObjArray *o(NULL);\r
281 const Char_t *tn=(name ? name : GetName());\r
282 if(!(o = (TObjArray*)gDirectory->Get(tn))){\r
283 AliWarning(Form("Missing histogram container %s.", tn));\r
284 return kFALSE;\r
285 }\r
286 fContainer = (TObjArray*)o->Clone(GetName());\r
287 gFile->Close();\r
288 return kTRUE;\r
289}\r
290\r
291//____________________________________________________________________\r
292void AliTRDinfoGen::UserExec(Option_t *){\r
293 //\r
294 // Run the Analysis\r
295 //\r
296\r
297 fTracksBarrel->Delete();\r
298 fTracksSA->Delete();\r
299 fTracksKink->Delete();\r
300 fV0List->Delete();\r
301 fEventInfo->Delete("");\r
302\r
303 fESDev = dynamic_cast<AliESDEvent*>(InputEvent());\r
304 if(!fESDev){\r
305 AliError("Failed retrieving ESD event");\r
306 return;\r
307 }\r
308 // WARNING\r
309 // This part may conflict with other detectors !!\r
4e38ea4c 310 if(!IsInitOCDB()){\r
94b94be0 311 AliCDBEntry* obj(NULL);\r
4e38ea4c 312 AliCDBManager* ocdb = AliCDBManager::Instance();\r
313 if(ocdb->IsDefaultStorageSet()){\r
314 AliInfo("OCDB :: initialized.");\r
94b94be0 315 } else {\r
4e38ea4c 316 AliInfo("OCDB :: initializing locally ...");\r
317 // prepare OCDB access\r
318 ocdb->SetDefaultStorage(fOCDB.Data());\r
319 ocdb->SetRun(fESDev->GetRunNumber());\r
320 // create geo manager\r
321 if(!(obj = ocdb->Get(AliCDBPath("GRP", "Geometry", "Data")))){\r
322 AliError("GEOMETRY failed initialization.");\r
323 } else {\r
324 AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());\r
325 AliGeomManager::GetNalignable("TRD");\r
326 AliGeomManager::ApplyAlignObjsFromCDB("TRD");\r
327 }\r
328 //init magnetic field\r
329 if(!TGeoGlobalMagField::Instance()->IsLocked() &&\r
330 !fESDev->InitMagneticField()){\r
331 AliError("MAGNETIC FIELD failed initialization.");\r
332 }\r
333 // set no of time bins\r
334 AliTRDtrackerV1::SetNTimeBins(AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS());\r
94b94be0 335 }\r
4e38ea4c 336 AliInfo(Form("OCDB : Loc[%s] Run[%d] TB[%d]", ocdb->GetDefaultStorage()->GetURI().Data(), ocdb->GetRun(), AliTRDtrackerV1::GetNTimeBins()));\r
337\r
338 // load misalignment\r
94b94be0 339 fgGeo = new AliTRDgeometry;\r
340 fgGeo->CreateClusterMatrixArray();\r
246fe7f7 341 // load reco param list from OCDB\r
99cd1793 342 AliInfo("Initializing TRD reco params ...");\r
94b94be0 343 fgReconstructor = new AliTRDReconstructor();\r
344 if(!(obj = ocdb->Get(AliCDBPath("TRD", "Calib", "RecoParam")))){\r
345 AliError("RECO PARAM failed initialization.");\r
94b94be0 346 } else {\r
347 obj->PrintMetaData();\r
246fe7f7 348 fRecos = (TObjArray*)obj->GetObject();\r
94b94be0 349 }\r
350 SetInitOCDB();\r
351 }\r
3ceb45ae 352\r
353 AliDebug(2, "+++++++++++++++++++++++++++++++++++++++");\r
354 AliDebug(2, Form("Analyse Ev[%d]", fESDev->GetEventNumberInFile()));\r
355\r
246fe7f7 356 // set reco param valid for this event\r
3ceb45ae 357 TH1 *h = (TH1I*)fContainer->At(kEvType);\r
246fe7f7 358 if(!fRecos){\r
359 fgReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());\r
360 h->Fill(0);\r
361 } else {\r
362 for(Int_t ireco(0); ireco<fRecos->GetEntriesFast(); ireco++){\r
363 AliTRDrecoParam *reco((AliTRDrecoParam*)fRecos->At(ireco));\r
364 Int_t es(reco->GetEventSpecie());\r
365 if(!(es&fESDev->GetEventSpecie())) continue;\r
366 fgReconstructor->SetRecoParam(reco);\r
3ceb45ae 367 if(AliLog::GetDebugLevel("PWG1/TRD", "AliTRDinfoGen")>5) reco->Dump();\r
246fe7f7 368 TString s;\r
369 if(es&AliRecoParam::kLowMult){ s="LowMult"; h->Fill(0);}\r
370 else if(es&AliRecoParam::kHighMult){ s="HighMult"; h->Fill(1);}\r
371 else if(es&AliRecoParam::kCosmic){ s="Cosmic"; h->Fill(2);}\r
372 else if(es&AliRecoParam::kCalib){ s="Calib"; h->Fill(3);}\r
373 else s="Unknown";\r
3ceb45ae 374 AliDebug(2, Form(" Using reco param \"%s\" for event %d.", s.Data(), fESDev->GetEventNumberInFile()));\r
246fe7f7 375 break;\r
376 }\r
377 }\r
94b94be0 378\r
379 // link MC if available\r
380 fMCev = MCEvent();\r
381\r
382 // event selection : trigger cut\r
383 if(UseLocalEvSelection() && fEvTrigger){ \r
384 Bool_t kTRIGGERED(kFALSE);\r
385 std::auto_ptr<TObjArray> trig(fEvTrigger->Tokenize(" "));\r
386 for(Int_t itrig=trig->GetEntriesFast(); itrig--;){\r
387 const Char_t *trigClass(((TObjString*)(*trig)[itrig])->GetName());\r
388 if(fESDev->IsTriggerClassFired(trigClass)) {\r
389 AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));\r
390 kTRIGGERED = kTRUE;\r
391 break; \r
392 }\r
393 }\r
394 if(!kTRIGGERED){ \r
395 AliDebug(2, Form("Reject Ev[%4d] Trigger", fESDev->GetEventNumberInFile()));\r
396 return;\r
397 }\r
398 // select only physical events\r
399 if(fESDev->GetEventType() != 7){ \r
400 AliDebug(2, Form("Reject Ev[%4d] EvType[%d]", fESDev->GetEventNumberInFile(), fESDev->GetEventType()));\r
401 return;\r
402 }\r
403 }\r
404\r
405 // if the required trigger is a collision trigger then apply event vertex cut\r
406 if(UseLocalEvSelection() && IsCollision()){\r
407 const AliESDVertex *vertex = fESDev->GetPrimaryVertex();\r
408 if(TMath::Abs(vertex->GetZv())<1.e-10 || \r
409 TMath::Abs(vertex->GetZv())>fgkEvVertexZ || \r
410 vertex->GetNContributors()<fgkEvVertexN) {\r
411 AliDebug(2, Form("Reject Ev[%4d] Vertex Zv[%f] Nv[%d]", fESDev->GetEventNumberInFile(), TMath::Abs(vertex->GetZv()), vertex->GetNContributors()));\r
412 return;\r
413 }\r
414 }\r
415\r
416 if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;\r
417\r
418 if(!fESDfriend){\r
419 AliError("Failed retrieving ESD friend event");\r
420 return;\r
421 }\r
422 if(HasMCdata() && !fMCev){\r
423 AliError("Failed retrieving MC event");\r
424 return;\r
425 }\r
426\r
427 new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));\r
428 // Determine centrality\r
429 // Author: Ionut Arsene <I.C.Arsene@gsi.de>\r
430 Int_t centralityBin = -1;\r
3ceb45ae 431 AliDebug(2, Form(" Beam Type: %s", fESDev->GetESDRun()->GetBeamType()));\r
94b94be0 432 if(TString(fESDev->GetESDRun()->GetBeamType()).Contains("Pb-Pb")){\r
433 centralityBin = 4;\r
434 const AliMultiplicity *mult = fESDev->GetMultiplicity();\r
435 Double_t zdcNeutronEnergy = fESDev->GetZDCN1Energy()+fESDev->GetZDCN2Energy();\r
436 Double_t itsNTracklets = mult->GetNumberOfTracklets();\r
437 Double_t centralitySlopes[6] = {0.0, 4.0, 8.0, 20.0, 50.0, 1000000.};\r
438 AliDebug(1, Form("zdcNeutronEnergy: %f, itsNTracklets: %f\n", zdcNeutronEnergy, itsNTracklets));\r
439 for(Int_t iCent=1; iCent<=5; ++iCent) {\r
440 if(zdcNeutronEnergy>centralitySlopes[iCent-1]*itsNTracklets && zdcNeutronEnergy<centralitySlopes[iCent]*itsNTracklets)\r
441 centralityBin=iCent - 1;\r
442 }\r
3ceb45ae 443 AliDebug(2, Form(" Centrality Class: %d", centralityBin));\r
94b94be0 444 }\r
445 fEventInfo->SetCentrality(centralityBin);\r
3ceb45ae 446 AliDebug(2, Form(" Bunch Crosses: %d", fESDev->GetBunchCrossNumber()));\r
94b94be0 447\r
448 Bool_t *trackMap(NULL);\r
449 AliStack * mStack(NULL);\r
450 Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();\r
451 if(HasMCdata()){\r
452 mStack = fMCev->Stack();\r
453 if(!mStack){\r
454 AliError("Failed retrieving MC Stack");\r
455 return;\r
456 }\r
457 trackMap = new Bool_t[nTracksMC];\r
458 memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);\r
459 }\r
460 \r
461 Double32_t dedx[100]; Int_t nSlices(0);\r
462 Int_t nTRDout(0), nTRDin(0), nTPC(0)\r
463 ,nclsTrklt\r
464 ,nBarrel(0), nSA(0), nKink(0)\r
465 ,nBarrelFriend(0), nSAFriend(0)\r
466 ,nBarrelMC(0), nSAMC(0), nKinkMC(0);\r
467 AliESDtrack *esdTrack = NULL;\r
468 AliESDfriendTrack *esdFriendTrack = NULL;\r
469 TObject *calObject = NULL;\r
470 AliTRDtrackV1 *track = NULL;\r
471 AliTRDseedV1 *tracklet = NULL;\r
472 AliTRDcluster *cl = NULL;\r
473\r
474\r
475 // LOOP 0 - over ESD v0s\r
476 Float_t bField(fESDev->GetMagneticField());\r
477 AliESDv0 *v0(NULL);\r
478 Int_t v0pid[AliPID::kSPECIES];\r
479 for(Int_t iv0(0); iv0<fESDev->GetNumberOfV0s(); iv0++){\r
480 if(!(v0 = fESDev->GetV0(iv0))) continue;\r
481 // register v0\r
482 if(fV0Cut) new(fV0Info) AliTRDv0Info(*fV0Cut);\r
483 else new(fV0Info) AliTRDv0Info();\r
484 fV0Info->SetMagField(bField);\r
485 fV0Info->SetV0tracks(fESDev->GetTrack(v0->GetPindex()), fESDev->GetTrack(v0->GetNindex()));\r
486 fV0Info->SetV0Info(v0);\r
487 fV0List->Add(new AliTRDv0Info(*fV0Info));// kFOUND=kFALSE;\r
488 }\r
489\r
490\r
491 // LOOP 1 - over ESD tracks\r
492 AliTRDv0Info *v0info=NULL;\r
493 for(Int_t itrk = 0; itrk < nTracksESD; itrk++){\r
494 new(fTrackInfo) AliTRDtrackInfo();\r
495 esdTrack = fESDev->GetTrack(itrk);\r
3ceb45ae 496 AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d] TOF-BC[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2), esdTrack->GetTOFBunchCrossing()));\r
94b94be0 497\r
498 if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;\r
499 if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;\r
500 if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;\r
501\r
502 // look at external track param\r
503 const AliExternalTrackParam *op = esdTrack->GetOuterParam();\r
504 Double_t xyz[3];\r
505 if(op){\r
506 op->GetXYZ(xyz);\r
507 op->Global2LocalPosition(xyz, op->GetAlpha());\r
508 AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));\r
509 }\r
510\r
511 // read MC info\r
512 Int_t fPdg = -1;\r
513 Int_t label = -1; UInt_t alab=UINT_MAX;\r
514 if(HasMCdata()){\r
515 label = esdTrack->GetLabel(); \r
516 alab = TMath::Abs(label);\r
517 // register the track\r
518 if(alab < UInt_t(nTracksMC)){ \r
519 trackMap[alab] = kTRUE; \r
520 } else { \r
521 AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));\r
522 continue; \r
523 }\r
524 AliMCParticle *mcParticle = NULL; \r
525 if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){\r
526 AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));\r
527 continue;\r
528 }\r
529 fPdg = mcParticle->Particle()->GetPdgCode();\r
530 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();\r
531 Int_t iref = 0; AliTrackReference *ref = NULL; \r
532 while(iref<nRefs){\r
533 ref = mcParticle->GetTrackReference(iref);\r
534 if(ref->LocalX() > fgkTPC) break;\r
535 iref++;\r
536 }\r
537\r
538 fTrackInfo->SetMC();\r
539 fTrackInfo->SetPDG(fPdg);\r
540 fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());\r
541 fTrackInfo->SetLabel(label);\r
542 Int_t jref = iref;//, kref = 0;\r
543 while(jref<nRefs){\r
544 ref = mcParticle->GetTrackReference(jref);\r
545 if(ref->LocalX() > fgkTRD) break;\r
546 AliDebug(4, Form(" trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));\r
547 fTrackInfo->AddTrackRef(ref);\r
548 jref++;\r
549 }\r
550 AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));\r
551 }\r
552\r
553 // copy some relevant info to TRD track info\r
554 fTrackInfo->SetStatus(esdTrack->GetStatus());\r
555 fTrackInfo->SetTrackId(esdTrack->GetID());\r
556 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);\r
557 fTrackInfo->SetESDpid(p);\r
558 fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());\r
559 if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();\r
560 memset(dedx, 0, 100*sizeof(Double32_t));\r
561 Int_t in(0);\r
562 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)\r
563 for(Int_t is=0; is<nSlices; is++) \r
564 dedx[in++]=esdTrack->GetTRDslice(il, is);\r
565 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);\r
566 fTrackInfo->SetSlices(in, dedx);\r
567 fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));\r
568 // some other Informations which we may wish to store in order to find problematic cases\r
569 fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));\r
570 fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));\r
3ceb45ae 571 fTrackInfo->SetTOFbc(static_cast<UShort_t>(esdTrack->GetTOFBunchCrossing()));\r
94b94be0 572 nclsTrklt = 0;\r
573 \r
574 // set V0pid info\r
575 for(Int_t iv(0); iv<fV0List->GetEntriesFast(); iv++){\r
576 if(!(v0info = (AliTRDv0Info*)fV0List->At(iv))) continue;\r
577 if(!v0info->GetV0Daughter(1) && !v0info->GetV0Daughter(-1)) continue;\r
578 if(!v0info->HasTrack(fTrackInfo)) continue;\r
579 memset(v0pid, 0, AliPID::kSPECIES*sizeof(Int_t));\r
580 fTrackInfo->SetV0();\r
581 for(Int_t is=AliPID::kSPECIES; is--;){v0pid[is] = v0info->GetPID(is, fTrackInfo);}\r
582 fTrackInfo->SetV0pid(v0pid);\r
583 fTrackInfo->SetV0();\r
584 //const AliTRDtrackInfo::AliESDinfo *ei = fTrackInfo->GetESDinfo();\r
585 break;\r
586 }\r
587\r
588 // read REC info\r
9b524890 589 esdFriendTrack = (fESDfriend->GetNumberOfTracks() > itrk) ? fESDfriend->GetTrack(itrk): NULL;\r
590\r
94b94be0 591 if(esdFriendTrack){\r
592 Int_t icalib = 0;\r
593 while((calObject = esdFriendTrack->GetCalibObject(icalib++))){\r
594 if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack\r
595 if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;\r
596 AliDebug(4, Form("TRD track OK"));\r
597 // Set the clusters to unused\r
598 for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){\r
599 if(!(tracklet = track->GetTracklet(ipl))) continue;\r
600 tracklet->ResetClusterIter();\r
601 while((cl = tracklet->NextCluster())) cl->Use(0);\r
602 }\r
603 fTrackInfo->SetTrack(track);\r
3ceb45ae 604 h = (TH1I*)fContainer->At(kBunchCross); h->Fill(esdTrack->GetTOFBunchCrossing());\r
94b94be0 605 break;\r
606 }\r
607 AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));\r
608 } else AliDebug(3, "No ESD friends");\r
609 if(op) fTrackInfo->SetOuterParam(op);\r
610\r
611 if(DebugLevel() >= 1){\r
612 AliTRDtrackInfo info(*fTrackInfo);\r
613 (*DebugStream()) << "trackInfo"\r
614 << "TrackInfo.=" << &info\r
615 << "\n";\r
616 info.Delete("");\r
617 }\r
618\r
619 ULong_t status(esdTrack->GetStatus());\r
620 if((status&AliESDtrack::kTPCout)){\r
621 if(!esdTrack->GetKinkIndex(0)){ // Barrel Track Selection\r
622 Bool_t selected(kTRUE);\r
623 if(UseLocalTrkSelection()){\r
624 if(esdTrack->Pt() < fgkPt){ \r
3ceb45ae 625 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] Pt[%5.2f]", itrk, fESDev->GetEventNumberInFile(), esdTrack->Pt()));\r
94b94be0 626 selected = kFALSE;\r
627 }\r
628 if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){\r
3ceb45ae 629 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] Eta[%5.2f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(esdTrack->Eta())));\r
94b94be0 630 selected = kFALSE;\r
631 }\r
632 if(selected && esdTrack->GetTPCNcls() < fgkNclTPC){ \r
3ceb45ae 633 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] NclTPC[%d]", itrk, fESDev->GetEventNumberInFile(), esdTrack->GetTPCNcls()));\r
94b94be0 634 selected = kFALSE;\r
635 }\r
636 Float_t par[2], cov[3];\r
637 esdTrack->GetImpactParameters(par, cov);\r
638 if(IsCollision()){ // cuts on DCA\r
639 if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){ \r
3ceb45ae 640 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] DCAxy[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[0])));\r
94b94be0 641 selected = kFALSE;\r
642 }\r
643 if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){ \r
3ceb45ae 644 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] DCAz[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[1])));\r
94b94be0 645 selected = kFALSE;\r
646 }\r
647 } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){;\r
3ceb45ae 648 AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] Primary", itrk, fESDev->GetEventNumberInFile()));\r
94b94be0 649 selected = kFALSE;\r
650 }\r
651 }\r
652 if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;\r
653 if(selected){ \r
654 fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));\r
655 nBarrel++;\r
656 if(fTrackInfo->GetTrack()) \r
657 nBarrelFriend++;\r
658 }\r
659 } else {\r
660 fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));\r
661 nKink++;\r
662 }\r
663 } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)){ \r
664 fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));\r
665 nSA++;\r
666 if(fTrackInfo->GetTrack()) \r
667 nSAFriend++;\r
668 }\r
669 fTrackInfo->Delete("");\r
670 }\r
671\r
672 // LOOP 2 - over MC tracks which are passing TRD where the track is not reconstructed\r
673 if(HasMCdata()){\r
674 AliDebug(10, "Output of the MC track map:");\r
675 for(Int_t itk = 0; itk < nTracksMC; itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));\r
676 \r
677 for(Int_t itk = 0; itk < nTracksMC; itk++){\r
678 if(trackMap[itk]) continue;\r
679 AliMCParticle *mcParticle = (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));\r
680 Int_t fPdg = mcParticle->Particle()->GetPdgCode();\r
681 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();\r
682 Int_t iref = 0; AliTrackReference *ref = NULL; \r
683 Int_t nRefsTRD = 0;\r
684 new(fTrackInfo) AliTRDtrackInfo();\r
685 fTrackInfo->SetMC();\r
686 fTrackInfo->SetPDG(fPdg);\r
687 while(iref<nRefs){ // count TRD TR\r
688 Bool_t kIN(kFALSE);\r
689 ref = mcParticle->GetTrackReference(iref);\r
690 if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTRD){\r
691 fTrackInfo->AddTrackRef(ref);\r
692 nRefsTRD++;kIN=kTRUE;\r
693 }\r
694 AliDebug(4, Form(" trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));\r
695 iref++;\r
696 }\r
697 if(!nRefsTRD){\r
698 // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the \r
699 // analysis job\r
700 fTrackInfo->Delete("");\r
701 continue;\r
702 }\r
703 fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());\r
704 fTrackInfo->SetLabel(itk);\r
705 if(DebugLevel() >= 1){\r
706 AliTRDtrackInfo info(*fTrackInfo);\r
707 (*DebugStream()) << "trackInfo"\r
708 << "TrackInfo.=" << &info\r
709 << "\n";\r
710 info.Delete("");\r
711 }\r
712 AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));\r
713 // check where the track starts\r
714 ref = mcParticle->GetTrackReference(0);\r
715 if(ref->LocalX() < fgkITS){ \r
716 fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));\r
717 nBarrelMC++;\r
718 } else if(ref->LocalX() < fgkTPC) {\r
719 fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));\r
720 nKinkMC++;\r
721 } else if(nRefsTRD>6){\r
722 fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));\r
723 nSAMC++;\r
724 }\r
725 fTrackInfo->Delete("");\r
726 }\r
727 delete[] trackMap;\r
728 }\r
729 AliDebug(1, Form(\r
730 "\nEv[%3d] Tracks: ESD[%d] MC[%d] V0[%d]\n"\r
731 " TPCout[%d] TRDin[%d] TRDout[%d]\n"\r
732 " Barrel[%3d+%3d=%3d] SA[%2d+%2d=%2d] Kink[%2d+%2d=%2d]"\r
733 ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC, fV0List->GetEntries()\r
734 , nTPC, nTRDin, nTRDout\r
735 ,nBarrel, nBarrelMC, fTracksBarrel->GetEntries()\r
736 ,nSA, nSAMC, fTracksSA->GetEntries()\r
737 ,nKink, nKinkMC, fTracksKink->GetEntries()\r
738 ));\r
739 // save track statistics\r
3ceb45ae 740 h = (TH1I*)fContainer->At(kStatTrk);\r
94b94be0 741 h->Fill(Float_t(kTracksESD), nTracksESD);\r
742 h->Fill(Float_t(kTracksMC), nTracksMC);\r
743 h->Fill(Float_t(kV0), fV0List->GetEntries());\r
744 h->Fill(Float_t(kTPC), nTPC);\r
745 h->Fill(Float_t(kTRDin), nTRDin);\r
746 h->Fill(Float_t(kTRDout), nTRDout);\r
747 h->Fill(Float_t(kBarrel), nBarrel);\r
748 h->Fill(Float_t(kBarrelMC), nBarrelMC);\r
749 h->Fill(Float_t(kSA), nSA);\r
750 h->Fill(Float_t(kSAMC), nSAMC);\r
751 h->Fill(Float_t(kKink), nKink);\r
752 h->Fill(Float_t(kKinkMC), nKinkMC);\r
753 h->Fill(Float_t(kBarrelFriend), nBarrelFriend);\r
754 h->Fill(Float_t(kSAFriend), nSAFriend);\r
755\r
756 PostData(AliTRDpwg1Helper::kTracksBarrel, fTracksBarrel);\r
757 PostData(AliTRDpwg1Helper::kTracksSA, fTracksSA);\r
758 PostData(AliTRDpwg1Helper::kTracksKink, fTracksKink);\r
759 PostData(AliTRDpwg1Helper::kEventInfo, fEventInfo);\r
760 PostData(AliTRDpwg1Helper::kV0List, fV0List);\r
761}\r
762\r
763//____________________________________________________________________\r
764void AliTRDinfoGen::SetLocalV0Selection(const AliTRDv0Info *v0)\r
765{\r
766// Set V0 cuts from outside\r
767\r
768 if(!fV0Cut) fV0Cut = new AliTRDv0Info(*v0);\r
769 else new(fV0Cut) AliTRDv0Info(*v0);\r
770}\r
771\r
772//____________________________________________________________________\r
773void AliTRDinfoGen::SetTrigger(const Char_t *trigger)\r
774{\r
775 if(!fEvTrigger) fEvTrigger = new TString(trigger);\r
776 else (*fEvTrigger) = trigger;\r
777}\r
778\r
779//____________________________________________________________________\r
780TTreeSRedirector* AliTRDinfoGen::DebugStream()\r
781{\r
782// Manage debug stream for task\r
783 if(!fDebugStream){\r
784 TDirectory *savedir = gDirectory;\r
785 fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");\r
786 savedir->cd();\r
787 }\r
788 return fDebugStream;\r
789}\r
790\r
791\r
792\r