]>
Commit | Line | Data |
---|---|---|
a65a7e70 | 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: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */ | |
17 | ||
18 | //////////////////////////////////////////////////////////////////////////// | |
19 | // | |
20 | // Tender wagon for TRD performance/calibration train | |
21 | // | |
22 | // In this wagon the information from | |
23 | // - ESD | |
24 | // - Friends [if available] | |
25 | // - MC [if available] | |
26 | // are grouped into AliTRDtrackInfo objects and fed to worker tasks | |
27 | // | |
28 | // Authors: | |
29 | // Markus Fasel <M.Fasel@gsi.de> | |
30 | // Alexandru Bercuci <A.Bercuci@gsi.de> | |
31 | // | |
32 | //////////////////////////////////////////////////////////////////////////// | |
33 | ||
34 | #include <TClonesArray.h> | |
35 | #include <TObjArray.h> | |
36 | #include <TObject.h> | |
37 | #include <TString.h> | |
38 | #include <TVectorT.h> | |
39 | #include <TH1.h> | |
40 | #include <TH2.h> | |
41 | #include <TCanvas.h> | |
42 | #include <TPad.h> | |
43 | #include <TFile.h> | |
44 | #include <TTree.h> | |
45 | #include <TROOT.h> | |
46 | #include <TChain.h> | |
47 | #include <TParticle.h> | |
48 | ||
49 | #include "AliLog.h" | |
50 | #include "AliAnalysisManager.h" | |
51 | #include "AliGeomManager.h" | |
52 | #include "AliCDBManager.h" | |
53 | #include "AliCDBStorage.h" | |
54 | #include "AliCDBEntry.h" | |
55 | #include "AliCDBPath.h" | |
56 | #include "AliESDEvent.h" | |
57 | #include "AliMCEvent.h" | |
58 | #include "AliESDInputHandler.h" | |
59 | #include "AliESDInputHandlerRP.h" | |
60 | #include "AliMCEventHandler.h" | |
61 | ||
62 | #include "AliESDfriend.h" | |
63 | #include "AliESDfriendTrack.h" | |
64 | #include "AliESDHeader.h" | |
65 | #include "AliESDRun.h" | |
66 | #include "AliESDtrack.h" | |
67 | #include "AliESDv0.h" | |
68 | #include "AliESDtrackCuts.h" | |
69 | #include "AliESDv0KineCuts.h" | |
70 | #include "AliMCParticle.h" | |
71 | #include "AliMultiplicity.h" | |
72 | #include "AliCentrality.h" | |
73 | #include "AliPID.h" | |
74 | #include "AliStack.h" | |
75 | #include "AliTrackReference.h" | |
76 | #include "TTreeStream.h" | |
77 | ||
78 | #include <cstdio> | |
79 | #include <climits> | |
80 | #include <cstring> | |
81 | #include <iostream> | |
82 | #include <memory> | |
83 | ||
84 | #include "AliTRDReconstructor.h" | |
85 | #include "AliTRDrecoParam.h" | |
86 | #include "AliTRDcalibDB.h" | |
87 | #include "AliTRDtrackerV1.h" | |
88 | #include "AliTRDpadPlane.h" | |
89 | #include "AliTRDgeometry.h" | |
90 | #include "AliTRDtrackV1.h" | |
91 | #include "AliTRDseedV1.h" | |
92 | #include "AliTRDcluster.h" | |
93 | #include "AliTRDinfoGen.h" | |
94 | #include "AliTRDpwgppHelper.h" | |
95 | #include "info/AliTRDtrackInfo.h" | |
96 | #include "info/AliTRDeventInfo.h" | |
97 | #include "info/AliTRDv0Info.h" | |
98 | #include "info/AliTRDchmbInfo.h" | |
99 | #include "info/AliTRDtriggerInfo.h" | |
100 | #include "info/AliTRDeventCuts.h" | |
101 | ||
102 | ClassImp(AliTRDinfoGen) | |
103 | ||
104 | const Float_t AliTRDinfoGen::fgkITS = 100.; // to be checked | |
105 | const Float_t AliTRDinfoGen::fgkTPC = 290.; | |
106 | const Float_t AliTRDinfoGen::fgkTRD = 365.; | |
107 | ||
108 | const Float_t AliTRDinfoGen::fgkTrkDCAxy = 1.; | |
109 | const Float_t AliTRDinfoGen::fgkTrkDCAz = 1.; | |
110 | const Int_t AliTRDinfoGen::fgkNclTPC = 70; | |
111 | const Float_t AliTRDinfoGen::fgkPt = 0.2; | |
112 | const Float_t AliTRDinfoGen::fgkEta = 0.9; | |
113 | AliTRDReconstructor* AliTRDinfoGen::fgReconstructor(NULL); | |
114 | AliTRDgeometry* AliTRDinfoGen::fgGeo(NULL); | |
115 | //____________________________________________________________________ | |
116 | AliTRDinfoGen::AliTRDinfoGen() | |
117 | :AliAnalysisTaskSE() | |
118 | ,fESDev(NULL) | |
119 | ,fMCev(NULL) | |
120 | ,fEventCut(NULL) | |
121 | ,fTrackCut(NULL) | |
122 | ,fV0Identifier(NULL) | |
123 | ,fV0Cut(NULL) | |
124 | ,fOCDB("local://$ALICE_ROOT/OCDB") | |
125 | ,fTrackInfo(NULL) | |
126 | ,fEventInfo(NULL) | |
127 | ,fV0Info(NULL) | |
128 | ,fTracksBarrel(NULL) | |
129 | ,fTracksITS(NULL) | |
130 | ,fTracksSA(NULL) | |
131 | ,fTracksKink(NULL) | |
132 | ,fV0List(NULL) | |
133 | ,fClusters(NULL) | |
134 | ,fContainer(NULL) | |
135 | ,fRecos(NULL) | |
136 | ,fDebugStream(NULL) | |
137 | { | |
138 | // | |
139 | // Default constructor | |
140 | // | |
141 | SetNameTitle("TRDinfoGen", "MC-REC TRD-track list generator"); | |
142 | } | |
143 | ||
144 | //____________________________________________________________________ | |
145 | AliTRDinfoGen::AliTRDinfoGen(char* name) | |
146 | :AliAnalysisTaskSE(name) | |
147 | ,fESDev(NULL) | |
148 | ,fMCev(NULL) | |
149 | ,fEventCut(NULL) | |
150 | ,fTrackCut(NULL) | |
151 | ,fV0Identifier(NULL) | |
152 | ,fV0Cut(NULL) | |
153 | ,fOCDB("local://$ALICE_ROOT/OCDB") | |
154 | ,fTrackInfo(NULL) | |
155 | ,fEventInfo(NULL) | |
156 | ,fV0Info(NULL) | |
157 | ,fTracksBarrel(NULL) | |
158 | ,fTracksITS(NULL) | |
159 | ,fTracksSA(NULL) | |
160 | ,fTracksKink(NULL) | |
161 | ,fV0List(NULL) | |
162 | ,fClusters(NULL) | |
163 | ,fContainer(NULL) | |
164 | ,fRecos(NULL) | |
165 | ,fDebugStream(NULL) | |
166 | { | |
167 | // | |
168 | // Default constructor | |
169 | // | |
170 | SetTitle("MC-REC TRD-track list generator"); | |
171 | DefineOutput(AliTRDpwgppHelper::kTracksBarrel, TObjArray::Class()); | |
172 | DefineOutput(AliTRDpwgppHelper::kTracksITS, TObjArray::Class()); | |
173 | DefineOutput(AliTRDpwgppHelper::kTracksSA, TObjArray::Class()); | |
174 | DefineOutput(AliTRDpwgppHelper::kTracksKink, TObjArray::Class()); | |
175 | DefineOutput(AliTRDpwgppHelper::kEventInfo, AliTRDeventInfo::Class()); | |
176 | DefineOutput(AliTRDpwgppHelper::kV0List, TObjArray::Class()); | |
177 | DefineOutput(AliTRDpwgppHelper::kClusters, TObjArray::Class()); | |
178 | DefineOutput(AliTRDpwgppHelper::kMonitor, TObjArray::Class()); // histogram list | |
179 | } | |
180 | ||
181 | //____________________________________________________________________ | |
182 | AliTRDinfoGen::~AliTRDinfoGen() | |
183 | { | |
184 | // Destructor | |
185 | if(fgGeo) delete fgGeo; | |
186 | if(fgReconstructor) delete fgReconstructor; | |
187 | if(fDebugStream) delete fDebugStream; | |
188 | if(fV0Cut) delete fV0Cut; | |
189 | if(fV0Identifier) delete fV0Identifier; | |
190 | if(fTrackCut) delete fTrackCut; | |
191 | if(fEventCut) delete fEventCut; | |
192 | if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL; | |
193 | if(fEventInfo) delete fEventInfo; fEventInfo = NULL; | |
194 | if(fV0Info) delete fV0Info; fV0Info = NULL; | |
195 | if(fTracksBarrel){ | |
196 | fTracksBarrel->Delete(); delete fTracksBarrel; | |
197 | fTracksBarrel = NULL; | |
198 | } | |
199 | if(fTracksITS){ | |
200 | fTracksITS->Delete(); delete fTracksITS; | |
201 | fTracksITS = NULL; | |
202 | } | |
203 | if(fTracksSA){ | |
204 | fTracksSA->Delete(); delete fTracksSA; | |
205 | fTracksSA = NULL; | |
206 | } | |
207 | if(fTracksKink){ | |
208 | fTracksKink->Delete(); delete fTracksKink; | |
209 | fTracksKink = NULL; | |
210 | } | |
211 | if(fV0List){ | |
212 | fV0List->Delete(); | |
213 | delete fV0List; | |
214 | fV0List = NULL; | |
215 | } | |
216 | if(fClusters){ | |
217 | fClusters->Delete(); delete fClusters; | |
218 | fClusters = NULL; | |
219 | } | |
220 | if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){ | |
221 | fContainer->Delete(); | |
222 | delete fContainer; | |
223 | fContainer = NULL; | |
224 | } | |
225 | } | |
226 | ||
227 | //____________________________________________________________________ | |
228 | Bool_t AliTRDinfoGen::GetRefFigure(Int_t) | |
229 | { | |
230 | // General graphs for PWGPP/TRD train | |
231 | if(!gPad){ | |
232 | AliWarning("Please provide a canvas to draw results."); | |
233 | return kFALSE; | |
234 | } | |
235 | fContainer->At(kStatTrk)->Draw("bar"); | |
236 | return kTRUE; | |
237 | } | |
238 | ||
239 | //____________________________________________________________________ | |
240 | void AliTRDinfoGen::UserCreateOutputObjects() | |
241 | { | |
242 | // | |
243 | // Create Output Containers (TObjectArray containing 1D histograms) | |
244 | // | |
245 | ||
246 | fTrackInfo = new AliTRDtrackInfo(); | |
247 | fEventInfo = new AliTRDeventInfo(); | |
248 | fV0Info = new AliTRDv0Info(); | |
249 | fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE); | |
250 | fTracksITS = new TObjArray(20); fTracksITS->SetOwner(kTRUE); | |
251 | fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE); | |
252 | fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE); | |
253 | fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE); | |
254 | fClusters = new TObjArray(AliTRDgeometry::kNdet); fClusters->SetOwner(kTRUE); | |
255 | ||
256 | // define general monitor | |
257 | fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE); | |
258 | TH1 *h=new TH1I("hStat", "Run statistics;Observable;Entries", Int_t(kNObjects), -0.5, Float_t(kNObjects)-0.5); | |
259 | TAxis *ax(h->GetXaxis()); | |
260 | ax->SetBinLabel(Int_t(kTracksESD) + 1, "ESD"); | |
261 | ax->SetBinLabel(Int_t(kTracksMC) + 1, "MC"); | |
262 | ax->SetBinLabel(Int_t(kV0) + 1, "V0"); | |
263 | ax->SetBinLabel(Int_t(kTPC) + 1, "TPC"); | |
264 | ax->SetBinLabel(Int_t(kITS) + 1, "ITS"); | |
265 | ax->SetBinLabel(Int_t(kTRDin) + 1, "TRDin"); | |
266 | ax->SetBinLabel(Int_t(kTRDout) + 1, "TRDout"); | |
267 | ax->SetBinLabel(Int_t(kBarrel) + 1, "Barrel"); | |
268 | ax->SetBinLabel(Int_t(kBarrelMC) + 1, "BarrelMC"); | |
269 | ax->SetBinLabel(Int_t(kSA) + 1, "SA"); | |
270 | ax->SetBinLabel(Int_t(kSAMC) + 1, "SAMC"); | |
271 | ax->SetBinLabel(Int_t(kKink) + 1, "Kink"); | |
272 | ax->SetBinLabel(Int_t(kKinkMC) + 1, "KinkMC"); | |
273 | ax->SetBinLabel(Int_t(kBarrelFriend) + 1, "BFriend"); | |
274 | ax->SetBinLabel(Int_t(kSAFriend) + 1, "SFriend"); | |
275 | fContainer->AddAt(h, kStatTrk); | |
276 | h=new TH1I("hEv", "Run statistics;Event Class;Entries", 4, -0.5, 3.5); | |
277 | ax = h->GetXaxis(); | |
278 | ax->SetBinLabel(1, "Low"); | |
279 | ax->SetBinLabel(2, "High"); | |
280 | ax->SetBinLabel(3, "Cosmic"); | |
281 | ax->SetBinLabel(4, "Calib"); | |
282 | fContainer->AddAt(h, kEvType); | |
283 | TH2I* h2=new TH2I("hBCtrack", "Track Statistics;Fill Bunch;TOF BC;Entries", 3500, -0.5, 3499.5, 31, -10.5, 20.5); | |
284 | fContainer->AddAt(h2, kBC); | |
285 | fContainer->AddAt(new AliTRDtriggerInfo(), kTrigger); | |
286 | TObjArray *chmb = new TObjArray(AliTRDgeometry::kNdet); | |
287 | chmb->SetName("Chambers Status"); chmb->SetOwner(kTRUE); | |
288 | fContainer->AddAt(chmb, kChmb); | |
289 | ||
290 | PostData(AliTRDpwgppHelper::kTracksBarrel, fTracksBarrel); | |
291 | PostData(AliTRDpwgppHelper::kTracksITS, fTracksITS); | |
292 | PostData(AliTRDpwgppHelper::kTracksSA, fTracksSA); | |
293 | PostData(AliTRDpwgppHelper::kTracksKink, fTracksKink); | |
294 | PostData(AliTRDpwgppHelper::kEventInfo, fEventInfo); | |
295 | PostData(AliTRDpwgppHelper::kV0List, fV0List); | |
296 | PostData(AliTRDpwgppHelper::kClusters, fClusters); | |
297 | PostData(AliTRDpwgppHelper::kMonitor, fContainer); | |
298 | } | |
299 | ||
300 | //____________________________________________________________________ | |
301 | Bool_t AliTRDinfoGen::Load(const Char_t *file, const Char_t *dir, const Char_t *name) | |
302 | { | |
303 | // Load data from performance file | |
304 | ||
305 | if(!TFile::Open(file)){ | |
306 | AliWarning(Form("Couldn't open file %s.", file)); | |
307 | return kFALSE; | |
308 | } | |
309 | if(!gFile->cd(dir)){ | |
310 | AliWarning(Form("Couldn't cd to %s in %s.", dir, file)); | |
311 | gFile->Close(); | |
312 | return kFALSE; | |
313 | } | |
314 | const Char_t *tn=(name ? name : GetName()); | |
315 | if(!(fContainer = (TObjArray*)gDirectory->Get(tn))){ | |
316 | AliWarning(Form("Missing histogram container %s.", tn)); | |
317 | gFile->Close(); | |
318 | return kFALSE; | |
319 | } | |
320 | gFile->Close(); | |
321 | return kTRUE; | |
322 | } | |
323 | ||
324 | //____________________________________________________________________ | |
325 | void AliTRDinfoGen::UserExec(Option_t *){ | |
326 | // | |
327 | // Run the Analysis | |
328 | // | |
329 | ||
330 | fTracksBarrel->Delete(); | |
331 | fTracksITS->Delete(); | |
332 | fTracksSA->Delete(); | |
333 | fTracksKink->Delete(); | |
334 | fV0List->Delete(); | |
335 | fClusters->Delete(); | |
336 | fEventInfo->Delete(""); | |
337 | ||
338 | fESDev = dynamic_cast<AliESDEvent*>(InputEvent()); | |
339 | if(!fESDev){ | |
340 | AliError("Failed retrieving ESD event"); | |
341 | return; | |
342 | } | |
343 | // WARNING | |
344 | // This part may conflict with other detectors !! | |
345 | if(!IsInitOCDB()){ | |
346 | AliCDBEntry* obj(NULL); | |
347 | AliCDBManager* ocdb = AliCDBManager::Instance(); | |
348 | if(ocdb->IsDefaultStorageSet()){ | |
349 | AliInfo("OCDB :: initialized externally."); | |
350 | } else { | |
351 | AliInfo("OCDB :: initializing locally ..."); | |
352 | // prepare OCDB access | |
353 | ocdb->SetDefaultStorage(fOCDB.Data()); | |
354 | //ocdb->SetSpecificStorage("TRD/Align/Data", "local:///home/niham/abercuci/local"); | |
355 | ocdb->SetRun(fESDev->GetRunNumber()); | |
356 | // create geo manager | |
357 | if(!(obj = ocdb->Get(AliCDBPath("GRP", "Geometry", "Data")))){ | |
358 | AliError("GEOMETRY failed initialization."); | |
359 | } else { | |
360 | AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject()); | |
361 | AliGeomManager::GetNalignable("TRD"); | |
362 | AliGeomManager::GetNalignable("TPC"); | |
363 | AliGeomManager::GetNalignable("ITS"); | |
364 | AliGeomManager::ApplyAlignObjsFromCDB("ITS TPC TRD"); | |
365 | } | |
366 | //init magnetic field | |
367 | if(!TGeoGlobalMagField::Instance()->IsLocked() && | |
368 | !fESDev->InitMagneticField()){ | |
369 | AliError("MAGNETIC FIELD failed initialization."); | |
370 | } | |
371 | // set no of time bins | |
372 | AliTRDtrackerV1::SetNTimeBins(AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS()); | |
373 | } | |
374 | AliInfo(Form("OCDB : Loc[%s] Run[%d] TB[%d]", ocdb->GetDefaultStorage()->GetURI().Data(), ocdb->GetRun(), AliTRDtrackerV1::GetNTimeBins())); | |
375 | ||
376 | // load misalignment | |
377 | fgGeo = new AliTRDgeometry; | |
378 | fgGeo->CreateClusterMatrixArray(); | |
379 | MakeChambers(); | |
380 | // load reco param list from OCDB | |
381 | AliInfo("Initializing TRD reco params ..."); | |
382 | fgReconstructor = new AliTRDReconstructor(); | |
383 | if(!(obj = ocdb->Get(AliCDBPath("TRD", "Calib", "RecoParam")))){ | |
384 | AliError("RECO PARAM failed initialization."); | |
385 | } else { | |
386 | obj->PrintMetaData(); | |
387 | fRecos = (TObjArray*)obj->GetObject(); | |
388 | } | |
389 | SetInitOCDB(); | |
390 | } | |
391 | ||
392 | AliDebug(2, "+++++++++++++++++++++++++++++++++++++++"); | |
393 | AliDebug(2, Form("Analyse Ev[%d]", fESDev->GetEventNumberInFile())); | |
394 | ||
395 | // set reco param valid for this event | |
396 | TH1 *h = (TH1I*)fContainer->At(kEvType); | |
397 | if(!fRecos){ | |
398 | fgReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam()); | |
399 | h->Fill(0); | |
400 | } else { | |
401 | for(Int_t ireco(0); ireco<fRecos->GetEntriesFast(); ireco++){ | |
402 | AliTRDrecoParam *reco((AliTRDrecoParam*)fRecos->At(ireco)); | |
403 | Int_t es(reco->GetEventSpecie()); | |
404 | if(!(es&fESDev->GetEventSpecie())) continue; | |
405 | fgReconstructor->SetRecoParam(reco); | |
406 | if(AliLog::GetDebugLevel("PWGPP/TRD", "AliTRDinfoGen")>5) reco->Dump(); | |
407 | TString s; | |
408 | if(es&AliRecoParam::kLowMult){ s="LowMult"; h->Fill(0);} | |
409 | else if(es&AliRecoParam::kHighMult){ s="HighMult"; h->Fill(1);} | |
410 | else if(es&AliRecoParam::kCosmic){ s="Cosmic"; h->Fill(2);} | |
411 | else if(es&AliRecoParam::kCalib){ s="Calib"; h->Fill(3);} | |
412 | else s="Unknown"; | |
413 | AliDebug(2, Form(" Using reco param \"%s\" for event %d.", s.Data(), fESDev->GetEventNumberInFile())); | |
414 | break; | |
415 | } | |
416 | } | |
417 | ||
418 | // link MC if available | |
419 | fMCev = MCEvent(); | |
420 | ||
421 | // trigger monitor | |
422 | AliTRDtriggerInfo *ti = (AliTRDtriggerInfo*)fContainer->At(kTrigger); | |
423 | TObjArray *evTriggers = fESDev->GetFiredTriggerClasses().Tokenize(" "); | |
424 | //printf("Ev[%03d] Triggers[%s]\n", fESDev->GetEventNumberInFile(), fESDev->GetFiredTriggerClasses().Data()); | |
425 | for(Int_t iet(evTriggers->GetEntriesFast()); iet--;) ti->Add((*evTriggers)[iet]->GetName()); | |
426 | evTriggers->Delete(); delete evTriggers; | |
427 | ||
428 | // event selection based on vertex cuts and trigger | |
429 | if(UseLocalEvSelection() && !fEventCut->IsSelected(fESDev, IsCollision())){ | |
430 | AliDebug(2, "Event failed selection on vertex and trigger"); | |
431 | return; | |
432 | } | |
433 | ||
434 | if(!fESDfriend){ | |
435 | AliError("Failed retrieving ESD friend event"); | |
436 | return; | |
437 | } | |
438 | if(HasMCdata() && !fMCev){ | |
439 | AliError("Failed retrieving MC event"); | |
440 | return; | |
441 | } | |
442 | ||
443 | new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun())); | |
444 | // Determine centrality | |
445 | // Author: Ionut Arsene <I.C.Arsene@gsi.de> | |
446 | TString beamtype = fESDev->GetESDRun()->GetBeamType(); | |
447 | if(beamtype.Contains("Pb-Pb") || beamtype.Contains("A-A")){ | |
448 | const AliMultiplicity *mult = fESDev->GetMultiplicity(); | |
449 | fEventInfo->SetMultiplicity(mult?mult->GetNumberOfTracklets():0); | |
450 | const AliCentrality *cent = fESDev->GetCentrality(); | |
451 | // centrality for different options V0 = "V0M", ITS = "TKL" etc | |
452 | fEventInfo->SetCentrality(cent?cent->GetCentralityPercentile("TKL"):-1.); | |
453 | AliDebug(2, Form(" Beam Type[%s] Mult[%d/%4d] Cent[%d/%5.2f]", | |
454 | beamtype.Data(), | |
455 | fEventInfo->GetMultiplicity(), mult?mult->GetNumberOfTracklets():0, | |
456 | fEventInfo->GetCentrality(), cent?cent->GetCentralityPercentile("TKL"):-1.)); | |
457 | } else { | |
458 | fEventInfo->SetMultiplicity(0); | |
459 | fEventInfo->SetCentrality(-1.); | |
460 | AliDebug(2, Form(" Beam Type[%s]", beamtype.Data())); | |
461 | } | |
462 | UShort_t evBC(fESDev->GetBunchCrossNumber()); | |
463 | ||
464 | // electron identifier from conversions | |
465 | if(!fV0Identifier) fV0Identifier = new AliESDv0KineCuts(); | |
466 | fV0Identifier->SetEvent(fESDev); | |
467 | ||
468 | Bool_t *trackMap(NULL); | |
469 | AliStack * mStack(NULL); | |
470 | Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks(); | |
471 | if(HasMCdata()){ | |
472 | mStack = fMCev->Stack(); | |
473 | if(!mStack){ | |
474 | AliError("Failed retrieving MC Stack"); | |
475 | return; | |
476 | } | |
477 | trackMap = new Bool_t[nTracksMC]; | |
478 | memset(trackMap, 0, sizeof(Bool_t) * nTracksMC); | |
479 | } | |
480 | ||
481 | Double32_t dedx[100]; Int_t nSlices(0); | |
482 | Int_t nTRDout(0), nTRDin(0), nTPC(0), nITS(0) | |
483 | // ,nclsTrklt | |
484 | ,nBarrel(0), nBarrelITS(0), nSA(0), nKink(0) | |
485 | ,nBarrelFriend(0), nBarrelITSFriend(0), nSAFriend(0) | |
486 | ,nBarrelMC(0), nSAMC(0), nKinkMC(0); | |
487 | AliESDtrack *esdTrack = NULL; | |
488 | AliESDfriendTrack *esdFriendTrack = NULL; | |
489 | TObject *calObject = NULL; | |
490 | AliTRDtrackV1 *track = NULL; | |
491 | AliTRDseedV1 *tracklet = NULL; | |
492 | AliTRDcluster *cl = NULL; | |
493 | ||
494 | ||
495 | // LOOP 0 - over ESD v0s | |
496 | Float_t bField(fESDev->GetMagneticField()); | |
497 | AliESDv0 *v0(NULL); | |
498 | Int_t v0pid[AliPID::kSPECIES]; | |
499 | for(Int_t iv0(0); iv0<fESDev->GetNumberOfV0s(); iv0++){ | |
500 | // Take only V0s from the On-the-fly v0 finder | |
501 | if(!(v0 = fESDev->GetV0(iv0))) continue; | |
502 | if(!v0->GetOnFlyStatus()) continue; | |
503 | // register v0 | |
504 | if(fV0Cut) new(fV0Info) AliTRDv0Info(*fV0Cut); | |
505 | else new(fV0Info) AliTRDv0Info(); | |
506 | fV0Info->SetMagField(bField); | |
507 | fV0Info->SetV0tracks(fESDev->GetTrack(v0->GetPindex()), fESDev->GetTrack(v0->GetNindex())); | |
508 | fV0Info->SetV0Info(v0); | |
509 | // tag conversion electrons | |
510 | Int_t pdgV0, pdgP, pdgN; | |
511 | if(fV0Identifier->ProcessV0(v0, pdgV0, pdgP, pdgN)) { | |
512 | switch(pdgV0){ | |
513 | case kGamma: fV0Info->SetDecay(AliTRDv0Info::kGamma); break; | |
514 | case kK0Short: fV0Info->SetDecay(AliTRDv0Info::kK0s); break; | |
515 | case kLambda0: fV0Info->SetDecay(AliTRDv0Info::kLambda); break; | |
516 | case kLambda0Bar: fV0Info->SetDecay(AliTRDv0Info::kAntiLambda); break; | |
517 | default: AliDebug(1, Form("V0[%+4d] -> +[%+4d] -[%+4d]. Decay not mapped.", pdgV0, pdgP, pdgN)); | |
518 | } | |
519 | } | |
520 | fV0List->Add(new AliTRDv0Info(*fV0Info));// kFOUND=kFALSE; | |
521 | } | |
522 | ||
523 | ||
524 | // LOOP 1 - over ESD tracks | |
525 | AliTRDv0Info *v0info=NULL; | |
526 | for(Int_t itrk = 0; itrk < nTracksESD; itrk++){ | |
527 | new(fTrackInfo) AliTRDtrackInfo(); | |
528 | esdTrack = fESDev->GetTrack(itrk); | |
529 | 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())); | |
530 | if(esdTrack->GetStatus()&AliESDtrack::kITSout) nITS++; | |
531 | if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++; | |
532 | if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++; | |
533 | if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++; | |
534 | ||
535 | /* Int_t ns(esdTrack->GetNumberOfTRDslices()); | |
536 | printf(" %3d ITS[%c] TPC[%c] TRDin[%c] TRDout[%c] TRDStop[%c] ns[%d]\n", itrk, | |
537 | (esdTrack->GetStatus()&AliESDtrack::kITSout)?'y':'n', | |
538 | (esdTrack->GetStatus()&AliESDtrack::kTPCout)?'y':'n', | |
539 | (esdTrack->GetStatus()&AliESDtrack::kTRDin)?'y':'n', | |
540 | (esdTrack->GetStatus()&AliESDtrack::kTRDout)?'y':'n', | |
541 | (esdTrack->GetStatus()&AliESDtrack::kTRDStop)?'y':'n', ns); | |
542 | if(ns){ | |
543 | for(Int_t ipl(0); ipl<AliTRDgeometry::kNlayer; ipl++){ | |
544 | Double_t sp, p(esdTrack->GetTRDmomentum(ipl, &sp)); | |
545 | printf(" [%d] p[%6.3f+-%6.3f] dEdx={", ipl, p, sp); | |
546 | for(Int_t is(0); is<8; is++) printf("%7.2f ", esdTrack->GetTRDslice(ipl, is)); printf("}\n"); | |
547 | } | |
548 | } | |
549 | */ | |
550 | // look at external track param | |
551 | const AliExternalTrackParam *op = esdTrack->GetOuterParam(); | |
552 | Double_t xyz[3]; | |
553 | if(op){ | |
554 | op->GetXYZ(xyz); | |
555 | op->Global2LocalPosition(xyz, op->GetAlpha()); | |
556 | AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0])); | |
557 | } | |
558 | ||
559 | // read MC info | |
560 | Int_t label = -1; UInt_t alab=UINT_MAX; | |
561 | if(HasMCdata()){ | |
562 | label = esdTrack->GetLabel(); | |
563 | alab = TMath::Abs(label); | |
564 | // register the track | |
565 | if(alab < UInt_t(nTracksMC)){ | |
566 | trackMap[alab] = kTRUE; | |
567 | } else { | |
568 | AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk)); | |
569 | continue; | |
570 | } | |
571 | AliMCParticle *mcParticle(NULL); | |
572 | if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){ | |
573 | AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk)); | |
574 | continue; | |
575 | } | |
576 | fTrackInfo->SetMC(); | |
577 | fTrackInfo->SetMCeta(mcParticle->Eta()); | |
578 | fTrackInfo->SetMCphi(mcParticle->Phi()); | |
579 | fTrackInfo->SetMCpt(mcParticle->Pt()); | |
580 | fTrackInfo->SetPDG(mcParticle->Particle()->GetPdgCode()); | |
581 | fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary()); | |
582 | fTrackInfo->SetLabel(label); | |
583 | fTrackInfo->SetTRDlabel(esdTrack->GetTRDLabel()); | |
584 | AliTrackReference *ref(NULL); | |
585 | for(Int_t iref(0); iref<mcParticle->GetNumberOfTrackReferences(); iref++){ | |
586 | if(!(ref = mcParticle->GetTrackReference(iref))) continue; | |
587 | if(ref->DetectorId() != AliTrackReference::kTRD) continue; | |
588 | AliDebug(4, Form(" TRD trackRef[%2d] @ r[%7.3f] [REC]", iref, ref->LocalX())); | |
589 | fTrackInfo->AddTrackRef(ref); | |
590 | } | |
591 | AliDebug(3, Form("Lab[%4d] pdg[%+4d] NtrackRefs[%d(%d)]", alab, mcParticle->Particle()->GetPdgCode(), fTrackInfo->GetNTrackRefs(), mcParticle->GetNumberOfTrackReferences())); | |
592 | } | |
593 | ||
594 | // copy some relevant info to TRD track info | |
595 | fTrackInfo->SetStatus(esdTrack->GetStatus()); | |
596 | fTrackInfo->SetTrackId(esdTrack->GetID()); | |
597 | Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p); | |
598 | fTrackInfo->SetESDpid(p); | |
599 | fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID()); | |
600 | fTrackInfo->SetESDeta(esdTrack->Eta()); | |
601 | Double_t loc[3]; | |
602 | if(esdTrack->GetXYZAt(298., fESDev->GetMagneticField(), loc)) fTrackInfo->SetESDphi(TMath::ATan2(loc[1], loc[0])); | |
603 | fTrackInfo->SetESDpt(esdTrack->Pt()); | |
604 | if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices(); | |
605 | memset(dedx, 0, 100*sizeof(Double32_t)); | |
606 | Int_t in(0); | |
607 | for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) | |
608 | for(Int_t is=0; is<nSlices; is++) | |
609 | dedx[in++]=esdTrack->GetTRDslice(il, is); | |
610 | for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il); | |
611 | fTrackInfo->SetSlices(in, dedx); | |
612 | fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2)); | |
613 | // some other Informations which we may wish to store in order to find problematic cases | |
614 | fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0)); | |
615 | fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1))); | |
616 | fTrackInfo->SetTPCdedx(esdTrack->GetTPCsignal()); | |
617 | Float_t tofTime = esdTrack->GetTOFsignal() - fESDev->GetT0TOF(0); | |
618 | fTrackInfo->SetTOFbeta(tofTime>0.?((esdTrack->GetIntegratedLength()/(tofTime*TMath::C()))*10e9):-999.); | |
619 | fTrackInfo->SetTOFbc(esdTrack->GetTOFBunchCrossing()==AliVTrack::kTOFBCNA?0:esdTrack->GetTOFBunchCrossing()); | |
620 | // nclsTrklt = 0; | |
621 | ||
622 | // set V0pid info | |
623 | //printf("%4d Looking for V0s...\n" , fTrackInfo->GetTrackId()); | |
624 | for(Int_t iv(0); iv<fV0List->GetEntriesFast(); iv++){ | |
625 | if(!(v0info = (AliTRDv0Info*)fV0List->At(iv))) continue; | |
626 | if(!v0info->GetV0Daughter(1) && !v0info->GetV0Daughter(-1)) continue; | |
627 | if(!v0info->HasTrack(fTrackInfo)) continue; | |
628 | //v0info->Print(); | |
629 | memset(v0pid, 0, AliPID::kSPECIES*sizeof(Int_t)); | |
630 | fTrackInfo->SetV0(); | |
631 | for(Int_t is=AliPID::kSPECIES; is--;) v0pid[is] = v0info->GetPID(is, fTrackInfo); fTrackInfo->SetV0pid(v0pid); | |
632 | if(v0info->IsDecay(AliTRDv0Info::kGamma)) fTrackInfo->SetElectron(); | |
633 | else if(v0info->IsDecay(AliTRDv0Info::kK0s)) fTrackInfo->SetPion(); | |
634 | else if(v0info->IsDecay(AliTRDv0Info::kLambda)) esdTrack->Charge()>0?fTrackInfo->SetProton():fTrackInfo->SetPion(); | |
635 | else if(v0info->IsDecay(AliTRDv0Info::kAntiLambda)) esdTrack->Charge()<0?fTrackInfo->SetProton():fTrackInfo->SetPion(); | |
636 | ||
637 | //TODO one track can be attached to more than one v0. Ideally one would need a list of v0 attached to the track info | |
638 | } | |
639 | ||
640 | // read track REC info | |
641 | if((esdFriendTrack = (fESDfriend->GetNumberOfTracks() > itrk) ? fESDfriend->GetTrack(itrk): NULL)) { | |
642 | fTrackInfo->SetTPCoutParam(esdFriendTrack->GetTPCOut()); | |
643 | fTrackInfo->SetITSoutParam(esdFriendTrack->GetITSOut()); | |
644 | const AliTrackPointArray *tps(NULL); | |
645 | if((tps=esdFriendTrack->GetTrackPointArray()) && HasTrackPoints()) fTrackInfo->SetTrackPointArray(tps); | |
646 | Int_t icalib = 0; | |
647 | while((calObject = esdFriendTrack->GetCalibObject(icalib++))){ | |
648 | if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack | |
649 | if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break; | |
650 | AliDebug(4, Form("TRD track OK")); | |
651 | // Set the clusters to unused | |
652 | for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){ | |
653 | if(!(tracklet = track->GetTracklet(ipl))) continue; | |
654 | tracklet->ResetClusterIter(); | |
655 | while((cl = tracklet->NextCluster())) cl->Use(0); | |
656 | } | |
657 | fTrackInfo->SetTrack(track); | |
658 | ((TH2I*)fContainer->At(kBC))->Fill(evBC, esdTrack->GetTOFBunchCrossing()); | |
659 | break; | |
660 | } | |
661 | AliDebug(3, Form("Ntracklets[%d]", fTrackInfo->GetNTracklets())); | |
662 | } else AliDebug(3, Form("No Friends for trk[%3d] Ntrk[%3d]", itrk, fESDfriend->GetNumberOfTracks())); | |
663 | if(op) fTrackInfo->SetOuterParam(op); | |
664 | ||
665 | if(DebugLevel() >= 1){ | |
666 | AliTRDtrackInfo info(*fTrackInfo); | |
667 | (*DebugStream()) << "trackInfo" | |
668 | << "TrackInfo.=" << &info | |
669 | << "\n"; | |
670 | info.Delete(""); | |
671 | } | |
672 | ||
673 | ULong_t status(esdTrack->GetStatus()); | |
674 | if((status&AliESDtrack::kTPCout)){ // TPC prolongation | |
675 | if(!esdTrack->GetKinkIndex(0)){ // Barrel Track Selection | |
676 | Bool_t selected(kTRUE); | |
677 | if(UseLocalTrkSelection()){ | |
678 | if(esdTrack->Pt() < fgkPt){ | |
679 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] Pt[%5.2f]", itrk, fESDev->GetEventNumberInFile(), esdTrack->Pt())); | |
680 | selected = kFALSE; | |
681 | } | |
682 | if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){ | |
683 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] Eta[%5.2f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(esdTrack->Eta()))); | |
684 | selected = kFALSE; | |
685 | } | |
686 | if(selected && esdTrack->GetTPCNcls() < fgkNclTPC){ | |
687 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] NclTPC[%d]", itrk, fESDev->GetEventNumberInFile(), esdTrack->GetTPCNcls())); | |
688 | selected = kFALSE; | |
689 | } | |
690 | if(selected && !(esdTrack->GetStatus()&AliESDtrack::kITSrefit)){ //SPD refit flag (Ionut) | |
691 | AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] ITSrefit", itrk, fESDev->GetEventNumberInFile())); | |
692 | selected = kFALSE; | |
693 | } | |
694 | UChar_t itsMap = esdTrack->GetITSClusterMap(); | |
695 | Bool_t firstSPDlayerHit = (itsMap & (UChar_t(1)<<0)), | |
696 | secondSPDlayerHit = (itsMap & (UChar_t(1)<<1)); | |
697 | if(selected && !(firstSPDlayerHit || secondSPDlayerHit)){ // request at least 1 SPD cluster (Ionut) | |
698 | AliDebug(3, Form("Reject Trk[%3d] Ev[%4d] no SPD cluster", itrk, fESDev->GetEventNumberInFile())); | |
699 | selected = kFALSE; | |
700 | } | |
701 | ||
702 | Float_t par[2], cov[3]; | |
703 | esdTrack->GetImpactParameters(par, cov); | |
704 | if(IsCollision()){ // cuts on DCA | |
705 | if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){ | |
706 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] DCAxy[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[0]))); | |
707 | selected = kFALSE; | |
708 | } | |
709 | if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){ | |
710 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] DCAz[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[1]))); | |
711 | selected = kFALSE; | |
712 | } | |
713 | } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){; | |
714 | AliDebug(3, Form("Reject TPC Trk[%3d] Ev[%4d] Primary", itrk, fESDev->GetEventNumberInFile())); | |
715 | selected = kFALSE; | |
716 | } | |
717 | } | |
718 | if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE; | |
719 | if(selected){ | |
720 | fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
721 | nBarrel++; | |
722 | if(fTrackInfo->GetTrack()) nBarrelFriend++; | |
723 | } | |
724 | } else { | |
725 | fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
726 | nKink++; | |
727 | } | |
728 | } else if((status&AliESDtrack::kITSout)) { // ITS prolongation | |
729 | Bool_t selected(kTRUE); | |
730 | if(UseLocalTrkSelection()){ | |
731 | if(esdTrack->Pt() < fgkPt){ | |
732 | AliDebug(3, Form("Reject ITS Trk[%3d] Ev[%4d] Pt[%5.2f]", itrk, fESDev->GetEventNumberInFile(), esdTrack->Pt())); | |
733 | selected = kFALSE; | |
734 | } | |
735 | if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){ | |
736 | AliDebug(3, Form("Reject ITS Trk[%3d] Ev[%4d] Eta[%5.2f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(esdTrack->Eta()))); | |
737 | selected = kFALSE; | |
738 | } | |
739 | Float_t par[2], cov[3]; | |
740 | esdTrack->GetImpactParameters(par, cov); | |
741 | if(IsCollision()){ // cuts on DCA | |
742 | if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){ | |
743 | AliDebug(3, Form("Reject ITS Trk[%3d] Ev[%4d] DCAxy[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[0]))); | |
744 | selected = kFALSE; | |
745 | } | |
746 | if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){ | |
747 | AliDebug(3, Form("Reject ITS Trk[%3d] Ev[%4d] DCAz[%f]", itrk, fESDev->GetEventNumberInFile(), TMath::Abs(par[1]))); | |
748 | selected = kFALSE; | |
749 | } | |
750 | } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){; | |
751 | AliDebug(3, Form("Reject ITS Trk[%3d] Ev[%4d] Primary", itrk, fESDev->GetEventNumberInFile())); | |
752 | selected = kFALSE; | |
753 | } | |
754 | } | |
755 | if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE; | |
756 | if(selected){ | |
757 | fTracksITS->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
758 | nBarrelITS++; | |
759 | if(fTrackInfo->GetTrack()) nBarrelITSFriend++; | |
760 | } | |
761 | } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)){ // TRD SA tracking | |
762 | fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
763 | nSA++; | |
764 | if(fTrackInfo->GetTrack()) nSAFriend++; | |
765 | } | |
766 | fTrackInfo->Delete(""); | |
767 | } | |
768 | // read clusters REC info | |
769 | TTree * treeR(NULL); AliESDInputHandlerRP *esdRPhandler(NULL); | |
770 | if(fInputHandler) esdRPhandler = dynamic_cast<AliESDInputHandlerRP*>(fInputHandler); | |
771 | if(esdRPhandler){ | |
772 | if((treeR = esdRPhandler->GetTreeR("TRD"))) { | |
773 | TObjArray *recPoints(NULL); | |
774 | if((treeR->GetBranch("TRDcluster"))){ | |
775 | treeR->SetBranchAddress("TRDcluster", &recPoints); | |
776 | for(Int_t idet(0); idet<treeR->GetEntries(); idet++){ | |
777 | treeR->GetEntry(idet); | |
778 | if(!recPoints->GetEntries()){ | |
779 | AliDebug(1, Form("Missing entry %d from TreeR", idet)); | |
780 | continue; | |
781 | } | |
782 | AliTRDcluster *c = (AliTRDcluster*)(*recPoints)[0]; | |
783 | if(!c){ | |
784 | AliDebug(1, Form("Missing first cluster in entry %d from TreeR", idet)); | |
785 | continue; | |
786 | } | |
787 | fClusters->AddAt(recPoints->Clone(Form("%03d", c->GetDetector())), c->GetDetector()); | |
788 | } | |
789 | } else AliDebug(3, "No TRDcluster branch"); | |
790 | } else AliDebug(3, "No RecPoints"); | |
791 | } else AliDebug(3, "No AliESDInputHandlerRP"); | |
792 | ||
793 | ||
794 | // LOOP 2 - over MC tracks which are passing TRD where the track is not reconstructed | |
795 | if(HasMCdata()){ | |
796 | AliDebug(10, "Output of the MC track map:"); | |
797 | for(Int_t itk = 0; itk < nTracksMC; itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE")); | |
798 | ||
799 | AliTrackReference *ref(NULL); | |
800 | for(Int_t itk = 0; itk < nTracksMC; itk++){ | |
801 | if(trackMap[itk]) continue; | |
802 | AliMCParticle *mcParticle = (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk)); | |
803 | ||
804 | Int_t nRefsTRD(0); | |
805 | for(Int_t iref(0); iref<mcParticle->GetNumberOfTrackReferences(); iref++){ // count TRD TR | |
806 | if(!(ref = mcParticle->GetTrackReference(iref))) continue; | |
807 | if(ref->DetectorId() != AliTrackReference::kTRD) continue; | |
808 | if(!nRefsTRD){ // build track info for this pure MC track | |
809 | new(fTrackInfo) AliTRDtrackInfo(); | |
810 | fTrackInfo->SetMC(); | |
811 | fTrackInfo->SetPDG(mcParticle->Particle()->GetPdgCode()); | |
812 | } | |
813 | AliDebug(4, Form(" TRD trackRef[%2d] @ r[%7.3f] [MC]", iref, ref->LocalX())); | |
814 | fTrackInfo->AddTrackRef(ref); | |
815 | } | |
816 | // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the | |
817 | // analysis job | |
818 | if(!nRefsTRD) continue; | |
819 | AliDebug(3, Form("Lab[%4d] pdg[%+4d] NtrackRefs[%d(%d)]", itk, mcParticle->Particle()->GetPdgCode(), fTrackInfo->GetNTrackRefs(), mcParticle->GetNumberOfTrackReferences())); | |
820 | fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary()); | |
821 | fTrackInfo->SetLabel(itk); | |
822 | if(DebugLevel() >= 1){ | |
823 | AliTRDtrackInfo info(*fTrackInfo); | |
824 | (*DebugStream()) << "trackInfo" | |
825 | << "TrackInfo.=" << &info | |
826 | << "\n"; | |
827 | info.Delete(""); | |
828 | } | |
829 | AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD)); | |
830 | // check where the track starts | |
831 | ref = mcParticle->GetTrackReference(0); | |
832 | if(ref->LocalX() < fgkITS){ | |
833 | fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
834 | //fTracksITS->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
835 | nBarrelMC++; | |
836 | } else if(ref->LocalX() < fgkTPC) { | |
837 | fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
838 | nKinkMC++; | |
839 | } else if(nRefsTRD>6){ | |
840 | fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo)); | |
841 | nSAMC++; | |
842 | } | |
843 | fTrackInfo->Delete(""); | |
844 | } | |
845 | delete[] trackMap; | |
846 | } | |
847 | AliDebug(1, Form( | |
848 | "\nEv[%3d] Tracks: ESD[%d] MC[%d] V0[%d]\n" | |
849 | " TPCout[%d] ITSout[%d] TRDin[%d] TRDout[%d]\n" | |
850 | " Barrel[%3d+%3d=%3d] ITS[%3d=%3d] SA[%2d+%2d=%2d] Kink[%2d+%2d=%2d]" | |
851 | ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC, fV0List->GetEntries() | |
852 | , nTPC, nITS, nTRDin, nTRDout | |
853 | ,nBarrel, nBarrelMC, fTracksBarrel->GetEntries() | |
854 | ,nBarrelITS, fTracksITS->GetEntries() | |
855 | ,nSA, nSAMC, fTracksSA->GetEntries() | |
856 | ,nKink, nKinkMC, fTracksKink->GetEntries() | |
857 | )); | |
858 | // save track statistics | |
859 | h = (TH1I*)fContainer->At(kStatTrk); | |
860 | h->Fill(Float_t(kTracksESD), nTracksESD); | |
861 | h->Fill(Float_t(kTracksMC), nTracksMC); | |
862 | h->Fill(Float_t(kV0), fV0List->GetEntries()); | |
863 | h->Fill(Float_t(kTPC), nTPC); | |
864 | h->Fill(Float_t(kITS), nITS); | |
865 | h->Fill(Float_t(kTRDin), nTRDin); | |
866 | h->Fill(Float_t(kTRDout), nTRDout); | |
867 | h->Fill(Float_t(kBarrel), nBarrel); | |
868 | h->Fill(Float_t(kBarrelMC), nBarrelMC); | |
869 | h->Fill(Float_t(kSA), nSA); | |
870 | h->Fill(Float_t(kSAMC), nSAMC); | |
871 | h->Fill(Float_t(kKink), nKink); | |
872 | h->Fill(Float_t(kKinkMC), nKinkMC); | |
873 | h->Fill(Float_t(kBarrelFriend), nBarrelFriend); | |
874 | h->Fill(Float_t(kSAFriend), nSAFriend); | |
875 | } | |
876 | ||
877 | ||
878 | //____________________________________________________________________ | |
879 | void AliTRDinfoGen::MakeChambers() | |
880 | { | |
881 | // Build chamber position and status | |
882 | if(!fContainer){ | |
883 | AliError("Missing container"); | |
884 | return; | |
885 | } | |
886 | AliTRDcalibDB *calib(AliTRDcalibDB::Instance()); | |
887 | if(!calib){ | |
888 | AliError("No access to calibration data"); | |
889 | return; | |
890 | } | |
891 | TObjArray *chmb = (TObjArray*)fContainer->At(kChmb); | |
892 | Int_t stat(0); | |
893 | Double_t alpha(0.), cs(-2.), sn(0.), pos[4]; | |
894 | for(Int_t isec(0); isec<AliTRDgeometry::kNsector; isec++){ | |
895 | alpha = (0.5+isec)*AliTRDgeometry::GetAlpha(); | |
896 | cs = TMath::Cos(alpha); | |
897 | sn = TMath::Sin(alpha); | |
898 | ||
899 | for(Int_t istk(0); istk<AliTRDgeometry::kNstack; istk++){ | |
900 | for(Int_t ilyr(0); ilyr<AliTRDgeometry::kNlayer; ilyr++){ | |
901 | Int_t idet(AliTRDgeometry::GetDetector(ilyr, istk, isec)); | |
902 | TGeoHMatrix *matrix(fgGeo->GetClusterMatrix(idet)); | |
903 | if(!matrix){ | |
904 | AliDebug(2, Form("Missing matrix for %03d [%02d_%d_%d]", idet, isec, istk, ilyr)); | |
905 | continue; | |
906 | } | |
907 | AliDebug(2, Form("Read info for %03d [%02d_%d_%d]", idet, isec, istk, ilyr)); | |
908 | AliTRDpadPlane *pp(fgGeo->GetPadPlane(ilyr, istk)); | |
909 | Double_t zm(0.5 * (pp->GetRow0() + pp->GetRowEnd())), | |
910 | loc0[] = {AliTRDgeometry::AnodePos(), pp->GetCol0(), zm-pp->GetRow0()}, | |
911 | loc1[] = {AliTRDgeometry::AnodePos(), pp->GetColEnd(), zm-pp->GetRowEnd()}, | |
912 | glb[3] = {1,1,1}; | |
913 | matrix->LocalToMaster(loc0, glb); | |
914 | Float_t phi = TMath::ATan2(glb[0]*sn + glb[1]*cs, glb[0]*cs - glb[1]*sn), | |
915 | tgl = glb[2]/glb[0]/TMath::Sqrt(1.+glb[1]*glb[1]/glb[0]/glb[0]), | |
916 | eta = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl)))); | |
917 | pos[0] = eta; pos[1] = phi; | |
918 | matrix->LocalToMaster(loc1, glb); | |
919 | phi = TMath::ATan2(glb[0]*sn + glb[1]*cs, glb[0]*cs - glb[1]*sn); | |
920 | tgl = glb[2]/glb[0]/TMath::Sqrt(1.+glb[1]*glb[1]/glb[0]/glb[0]); | |
921 | eta = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl)))); | |
922 | pos[2] = eta; pos[3] = phi; | |
923 | stat = 0; | |
924 | if(calib->IsChamberGood(idet)){ | |
925 | if(calib->IsHalfChamberNoData(idet, 0)) stat += 2; | |
926 | if(calib->IsHalfChamberNoData(idet, 1)) stat += 3; | |
927 | } else stat = 1; | |
928 | chmb->Add(new AliTRDchmbInfo(idet, stat, pos)); | |
929 | } | |
930 | } | |
931 | } | |
932 | } | |
933 | ||
934 | //____________________________________________________________________ | |
935 | void AliTRDinfoGen::MakeSummary() | |
936 | { | |
937 | // Build summary plots | |
938 | if(!fContainer){ | |
939 | AliError("Missing results"); | |
940 | return; | |
941 | } | |
942 | TH1 *h1(NULL); TVirtualPad *p(NULL); TCanvas *cOut(NULL); | |
943 | ||
944 | const Int_t nx(1024), ny(1024); | |
945 | cOut = new TCanvas("infoGenSummary", "Run Statistics", nx, ny); | |
946 | cOut->Divide(2,2, 1.e-5, 1.e-5); | |
947 | //========= | |
948 | p=cOut->cd(1);p->SetRightMargin(0.025);p->SetTopMargin(0.01);p->SetLogy(); | |
949 | h1 = (TH1*)fContainer->At(kStatTrk); | |
950 | h1->SetBarOffset(0.06); h1->SetBarWidth(0.88); h1->SetFillColor(kGreen); h1->SetFillStyle(3001); | |
951 | h1->Draw("bar1"); | |
952 | //========= | |
953 | p=cOut->cd(2);p->SetRightMargin(0.025);p->SetTopMargin(0.01);p->SetLogy(); | |
954 | h1 = (TH1*)fContainer->At(kEvType); | |
955 | h1->SetBarOffset(0.04); h1->SetBarWidth(0.92);h1->SetFillColor(kGreen); h1->SetFillStyle(3001); | |
956 | h1->Draw("bar1"); | |
957 | //========= | |
958 | p=cOut->cd(3);p->SetRightMargin(0.025);p->SetTopMargin(0.01);p->SetLogz(); | |
959 | TH2 *h2((TH2*)fContainer->At(kBC)); | |
960 | h1 = h2->ProjectionX("_px"); | |
961 | Int_t n(0); Int_t bins[3500]; | |
962 | for(Int_t ib(1); ib<=3500; ib++){ | |
963 | if(h1->GetBinContent(ib) < 1) continue; | |
964 | bins[n++] = ib; | |
965 | } | |
966 | delete h1; | |
967 | ||
968 | TAxis *ay(h2->GetYaxis()); | |
969 | TH2 *hs = new TH2I("hBC_Summary", Form("%s;%s;%s", h2->GetTitle(), h2->GetXaxis()->GetTitle(), h2->GetYaxis()->GetTitle()), | |
970 | n, -0.5, n-0.5, ay->GetNbins(), ay->GetXmin(), ay->GetXmax()); | |
971 | hs->SetLineColor(kBlack);hs->SetLineWidth(1); | |
972 | hs->SetMarkerColor(kRed); | |
973 | TAxis *ax(hs->GetXaxis()); ax->CenterTitle(); ax->SetTitleOffset(1.4); | |
974 | for(Int_t ib(0); ib<n; ib++){ | |
975 | ax->SetBinLabel(ib+1, Form("%d", bins[ib])); | |
976 | for(Int_t iy(1); iy<=ay->GetNbins(); iy++){ | |
977 | hs->SetBinContent(ib+1, iy, h2->GetBinContent(bins[ib], iy)); | |
978 | } | |
979 | } | |
980 | hs->Draw("textbox"); | |
981 | ||
982 | //========= | |
983 | p=cOut->cd(4); p->SetRightMargin(0.0215);p->SetLeftMargin(0.414);//p->SetLogz(); | |
984 | TObject *o = fContainer->At(kTrigger); | |
985 | if(o){ | |
986 | if((h1 = dynamic_cast<TH1I*>(o))) { | |
987 | h1->GetXaxis()->SetTitleOffset(6.5); h1->GetXaxis()->CenterTitle(); | |
988 | h1->SetFillStyle(3001);h1->SetFillColor(kGreen); | |
989 | h1->SetBarWidth(0.8);h1->SetBarOffset(0.1); | |
990 | ((TH1I*)o)->Draw("hbar2"); | |
991 | } else o->Draw(); | |
992 | } | |
993 | cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
994 | } | |
995 | ||
996 | //____________________________________________________________________ | |
997 | void AliTRDinfoGen::SetLocalEvSelection(const AliTRDeventCuts &ec) | |
998 | { | |
999 | // Set event cuts from outside | |
1000 | if(!fEventCut) fEventCut = new AliTRDeventCuts(ec); | |
1001 | else new(fEventCut) AliTRDeventCuts(ec); | |
1002 | fEventCut->Print(); | |
1003 | } | |
1004 | ||
1005 | //____________________________________________________________________ | |
1006 | void AliTRDinfoGen::SetLocalV0Selection(const AliTRDv0Info &v0) | |
1007 | { | |
1008 | // Set V0 cuts from outside | |
1009 | ||
1010 | if(!fV0Cut) fV0Cut = new AliTRDv0Info(v0); | |
1011 | else new(fV0Cut) AliTRDv0Info(v0); | |
1012 | fV0Cut->Print(); | |
1013 | } | |
1014 | ||
1015 | //____________________________________________________________________ | |
1016 | TTreeSRedirector* AliTRDinfoGen::DebugStream() | |
1017 | { | |
1018 | // Manage debug stream for task | |
1019 | if(!fDebugStream){ | |
1020 | TDirectory *savedir = gDirectory; | |
1021 | fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root"); | |
1022 | savedir->cd(); | |
1023 | } | |
1024 | return fDebugStream; | |
1025 | } | |
1026 | ||
1027 | //____________________________________________________________________ | |
1028 | void AliTRDinfoGen::Terminate(Option_t* /*option*/) | |
1029 | { | |
1030 | // Process run information | |
1031 | AliInfo(""); | |
1032 | if(!(fContainer = dynamic_cast<TObjArray *>(GetOutputData(AliTRDpwgppHelper::kMonitor)))) return; | |
1033 | AliInfo(Form("fContainer(%p)", (void*)fContainer)); | |
1034 | ||
1035 | AliTRDtriggerInfo* ti(NULL); | |
1036 | if(UseLocalEvSelection()){ | |
1037 | if(!(ti = (AliTRDtriggerInfo*)fContainer->At(kTrigger))) return; | |
1038 | for(Int_t ix(0); ix<ti->GetNTriggers(); ix++){ | |
1039 | if(fEventCut->CheckTrigger(ti->GetTrigger(ix))) ti->SetSelectTrigger(ix); | |
1040 | //ax->SetBinLabel(ix, Form("#color[2]{%s}", ax->GetBinLabel(ix))); | |
1041 | } | |
1042 | } | |
1043 | } | |
1044 |