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