Updated PID information in D+ ntuple (Giacomo)
[u/mrichter/AliRoot.git] / PWGPP / TRD / AliTRDrecoTask.cxx
CommitLineData
94b94be0 1///////////////////////////////////////////////////////\r
2//\r
3// Basic class for Performance/Calibration TRD tasks\r
4// \r
5// It performs generic tasks like :\r
6// - data file manegment\r
7// - reference container management\r
8// - debug container management\r
9// - interaction with AliAnalysisManager\r
10// - Plot functor loop\r
11//\r
12// Author: Alexandru Bercuci <A.Bercuci@gsi.de>, 10/09/2008\r
13//\r
14/////////////////////////////////////////////////////////\r
15\r
16#include "TClass.h"\r
17#include "TMethod.h"\r
18#include "TMethodCall.h"\r
19#include "TMethodArg.h"\r
20#include "TFile.h"\r
21#include "TChain.h"\r
22#include "TList.h"\r
23#include "TMap.h"\r
24#include "TH1.h"\r
25#include "TF1.h"\r
26#include "TObjArray.h"\r
27#include "TDirectory.h"\r
28#include "TTreeStream.h"\r
82e6e5dc 29#include "TBox.h"\r
5591f3d6 30#include "TLatex.h"\r
82e6e5dc 31#include "TVectorT.h"\r
94b94be0 32\r
33#include "AliLog.h"\r
34#include "AliAnalysisTask.h"\r
100f4577 35#include "AliAnalysisManager.h"\r
3ceb45ae 36#include "AliExternalTrackParam.h"\r
94b94be0 37\r
38#include "info/AliTRDeventInfo.h"\r
39#include "AliTRDrecoTask.h"\r
3ceb45ae 40#include "AliTRDtrackV1.h"\r
41#include "AliTRDpidUtil.h"\r
94b94be0 42\r
43ClassImp(AliTRDrecoTask)\r
44\r
45TList* AliTRDrecoTask::fgTrendPoint(NULL);\r
46TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
47//_______________________________________________________\r
48AliTRDrecoTask::AliTRDrecoTask()\r
49 : AliAnalysisTaskSE()\r
50 ,fNRefFigures(0)\r
82e6e5dc 51 ,fDets(NULL)\r
94b94be0 52 ,fContainer(NULL)\r
53 ,fEvent(NULL)\r
54 ,fTracks(NULL)\r
f073d500 55 ,fClusters(NULL)\r
56 ,fkClusters(NULL)\r
94b94be0 57 ,fkTrack(NULL)\r
58 ,fkMC(NULL)\r
59 ,fkESD(NULL)\r
3ceb45ae 60 ,fSpecies(-6)\r
61 ,fPt(-1.)\r
62 ,fPhi(0.)\r
63 ,fEta(0.)\r
94b94be0 64 ,fPlotFuncList(NULL)\r
f073d500 65 ,fDetFuncList(NULL)\r
94b94be0 66 ,fRunTerminate(kFALSE)\r
67{\r
68// Default constructor\r
69 snprintf(fNameId, 10, "no name");\r
70}\r
71\r
72//_______________________________________________________\r
73AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
74 : AliAnalysisTaskSE(name)\r
75 ,fNRefFigures(0)\r
82e6e5dc 76 ,fDets(NULL)\r
94b94be0 77 ,fContainer(NULL)\r
78 ,fEvent(NULL)\r
79 ,fTracks(NULL)\r
f073d500 80 ,fClusters(NULL)\r
81 ,fkClusters(NULL)\r
94b94be0 82 ,fkTrack(NULL)\r
83 ,fkMC(NULL)\r
84 ,fkESD(NULL)\r
3ceb45ae 85 ,fSpecies(-6)\r
86 ,fPt(-1.)\r
87 ,fPhi(0.)\r
88 ,fEta(0.)\r
94b94be0 89 ,fPlotFuncList(NULL)\r
f073d500 90 ,fDetFuncList(NULL)\r
94b94be0 91 ,fRunTerminate(kFALSE)\r
92{\r
93// Constructor for all derived performance tasks\r
94\r
95 SetTitle(title);\r
96 snprintf(fNameId, 10, "no name");\r
97 DefineInput (1, TObjArray::Class()); // track list\r
98 DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
f073d500 99 DefineInput (3, TObjArray::Class()); // cluster list object\r
94b94be0 100 DefineOutput(1, TObjArray::Class()); // histogram list\r
101}\r
102\r
103//_______________________________________________________\r
104AliTRDrecoTask::~AliTRDrecoTask() \r
105{\r
106\r
107 // Generic task destructor\r
108\r
109 AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
110 if(fgDebugStream){ \r
111 delete fgDebugStream;\r
112 fgDebugStream = NULL;\r
113 }\r
114\r
115 if(fPlotFuncList){\r
116 fPlotFuncList->Delete();\r
117 delete fPlotFuncList;\r
118 fPlotFuncList = NULL;\r
119 }\r
f073d500 120 if(fDetFuncList){\r
121 fDetFuncList->Delete();\r
122 delete fDetFuncList;\r
123 fDetFuncList = NULL;\r
124 }\r
94b94be0 125 \r
82e6e5dc 126 if(fDets){\r
127 if(fDets->IsOwner()) fDets->Delete();\r
128 delete fDets;\r
129 fDets = NULL;\r
130 }\r
131\r
f0857a6a 132 if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
94b94be0 133 if(fContainer->IsOwner()) fContainer->Delete();\r
134 delete fContainer;\r
135 fContainer = NULL;\r
136 }\r
137\r
138 if(fgTrendPoint){\r
139 TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
140 fgTrendPoint->Write();\r
141 delete fgTrendPoint;\r
142 fgTrendPoint=NULL;\r
143 gFile->Close();\r
144 }\r
145}\r
146\r
147//_______________________________________________________\r
148Int_t AliTRDrecoTask::GetNRefFigures() const \r
149{ \r
150 if(!fNRefFigures) AliWarning("No reference plots available.");\r
151 return fNRefFigures; \r
152} \r
153\r
154//_______________________________________________________\r
155void AliTRDrecoTask::UserCreateOutputObjects()\r
156{\r
157 if(!HasFunctorList()) InitFunctorList();\r
158 fContainer = Histos();\r
159 PostData(1, fContainer);\r
160}\r
161\r
162//_______________________________________________________\r
163void AliTRDrecoTask::UserExec(Option_t *)\r
164{\r
165// Loop over Plot functors published by particular tasks\r
166\r
f073d500 167 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));\r
168 fEvent = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
169 fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
94b94be0 170\r
171 if(!fPlotFuncList){\r
f073d500 172 AliWarning("No track functor list defined for the task");\r
94b94be0 173 return;\r
174 }\r
175 if(!fTracks) return;\r
176 if(!fTracks->GetEntriesFast()) return;\r
177 else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
178\r
178d284a 179 Int_t itrk(-1);\r
180 AliTRDtrackInfo *trackInfo(NULL);\r
94b94be0 181 TIter plotIter(fPlotFuncList);\r
182 TObjArrayIter trackIter(fTracks);\r
183 while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
178d284a 184 itrk++; fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
94b94be0 185 fkMC = trackInfo->GetMCinfo();\r
186 fkESD = trackInfo->GetESDinfo();\r
178d284a 187 if((fkTrack = trackInfo->GetTrack())){\r
188 // cache properties of the track at TRD entrance\r
189 // check input track status\r
190 AliExternalTrackParam *tin(NULL);\r
191 if(!(tin = fkTrack->GetTrackIn())) AliDebug(2, Form("Missing TRD track[%d] :: entry point.", itrk));\r
192 else {\r
193 fPt = tin->Pt();\r
194 fEta = tin->Eta();\r
195 Double_t xyz[3];\r
196 if(!tin->GetXYZ(xyz)) AliDebug(2, Form("Failed TRD track[%d] :: global track postion", itrk));\r
197 else fPhi = TMath::ATan2(xyz[1], xyz[0]);\r
198 fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
199 }\r
200 } else AliDebug(2, Form("Missing TRD track[%d].", itrk));\r
201\r
202 TMethodCall *plot(NULL);\r
94b94be0 203 plotIter.Reset();\r
178d284a 204 while((plot=dynamic_cast<TMethodCall*>(plotIter()))) plot->Execute(this);\r
94b94be0 205 }\r
f073d500 206 if(!fClusters) return;\r
207 if(!fDetFuncList){\r
178d284a 208 AliDebug(1, "No detector functor list defined for task");\r
f073d500 209 return;\r
210 }\r
211 TIter detIter(fDetFuncList);\r
212 for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
213 if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
214 TMethodCall *det(NULL);\r
215 detIter.Reset();\r
178d284a 216 while((det=dynamic_cast<TMethodCall*>(detIter()))) det->Execute(this);\r
f073d500 217 }\r
94b94be0 218}\r
219\r
220//_______________________________________________________\r
221Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
222{\r
223 AliWarning("Retrieving reference figures not implemented.");\r
224 return kFALSE;\r
225}\r
226\r
227//_______________________________________________________\r
228Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val)\r
229{\r
230// Generic publisher for trend values\r
231\r
232 if(!fgTrendPoint){\r
233 fgTrendPoint = new TList();\r
234 fgTrendPoint->SetOwner();\r
235 }\r
236 fgTrendPoint->AddLast(new TNamed(Form("%s_%s", GetName(), name), Form("%f", val)));\r
237 return kTRUE;\r
238}\r
239\r
240//_______________________________________________________\r
241void AliTRDrecoTask::InitFunctorList()\r
242{\r
243// Initialize list of functors\r
244\r
245 TClass *c = this->IsA();\r
3ed01fbe 246 if(fPlotFuncList) fPlotFuncList->Clear();\r
f073d500 247 if(fDetFuncList) fDetFuncList->Clear();\r
94b94be0 248\r
f073d500 249 TMethod *m(NULL);\r
94b94be0 250 TIter methIter(c->GetListOfMethods());\r
251 while((m=dynamic_cast<TMethod*>(methIter()))){\r
252 TString name(m->GetName());\r
f073d500 253 if(name.BeginsWith("Plot")){\r
254 if(!fPlotFuncList) fPlotFuncList = new TList();\r
255 fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
256 } else if(name.BeginsWith("Det")){\r
257 if(!fDetFuncList) fDetFuncList = new TList();\r
258 fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
259 }\r
94b94be0 260 }\r
261}\r
262\r
263//_______________________________________________________\r
264Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
265{\r
266// Generic container loader\r
267\r
268 if(!TFile::Open(file)){\r
269 AliWarning(Form("Couldn't open file %s.", file));\r
270 return kFALSE;\r
271 }\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 TObjArray *o = NULL;\r
277 if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
278 AliWarning("Missing histogram container.");\r
279 return kFALSE;\r
280 }\r
281 fContainer = (TObjArray*)o->Clone(GetName());\r
282 gFile->Close();\r
283 return kTRUE;\r
284}\r
285\r
286//________________________________________________________\r
82e6e5dc 287Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
288{\r
289// Load detector map.\r
290\r
291 if(!TFile::Open(file)){\r
292 AliWarning(Form("Couldn't open file %s.", file));\r
293 return kFALSE;\r
294 }\r
295 if(!gFile->cd(dir)){\r
296 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
297 return kFALSE;\r
298 }\r
299 TObjArray *info = NULL;\r
300 if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
301 AliWarning("Missing TRDinfoGen container.");\r
302 return kFALSE;\r
303 }\r
304 TObjArray *dets = (TObjArray*)info->FindObject("Chambers");\r
305 if(!dets){\r
306 AliWarning("Missing detector map from TRDinfoGen results.");\r
307 info->ls();\r
308 return kFALSE;\r
309 }\r
310 fDets = (TObjArray*)dets->Clone("Chambers");\r
311 gFile->Close();\r
312 return kTRUE;\r
313}\r
314\r
315\r
316//________________________________________________________\r
94b94be0 317Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
318 //\r
319 // Store the output graphs in a ROOT file\r
320 // Input TObject array will not be written as Key to the file,\r
321 // only content itself\r
322 //\r
323\r
324 TDirectory *cwd = gDirectory;\r
325 if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
326\r
327 TIterator *iter = results->MakeIterator();\r
328 TObject *inObject = NULL, *outObject = NULL;\r
329 while((inObject = iter->Next())){\r
330 outObject = inObject->Clone();\r
331 outObject->Write(NULL, TObject::kSingleKey);\r
332 }\r
333 delete iter;\r
334 gFile->Close(); delete gFile;\r
335 cwd->cd(); \r
336 return kTRUE;\r
337}\r
338\r
339//_______________________________________________________\r
340Bool_t AliTRDrecoTask::PostProcess()\r
341{\r
342// To be implemented by particular tasks\r
343\r
344 AliWarning("Post processing of reference histograms not implemented.");\r
345 return kTRUE;\r
346}\r
347\r
348//_______________________________________________________\r
5591f3d6 349void AliTRDrecoTask::MakeDetectorPlot(Int_t ly, const Option_t *opt)\r
82e6e5dc 350{\r
351// Draw chamber boundaries in eta/phi plots with misalignments\r
352// based on info collected by AliTRDinfoGen\r
353\r
354 if(!fDets){\r
355 AliWarning("Detector map and status not available.");\r
356 return;\r
357 }\r
5591f3d6 358 Float_t xmin(0.), xmax(0.);\r
82e6e5dc 359 TBox *gdet = new TBox();\r
360 gdet->SetLineWidth(kBlack);gdet->SetFillColor(kBlack);\r
361 Int_t style[] = {0, 3003};\r
362 for(Int_t idet(0); idet<540; idet++){\r
363 if(idet%6 != ly) continue;\r
364 TVectorF *det((TVectorF*)fDets->At(idet));\r
365 if(!det) continue;\r
82e6e5dc 366 Int_t iopt = Int_t((*det)[4]);\r
b9058d72 367 if(strcmp(opt, "eta")==0){\r
368 xmin=(*det)[0]; xmax=(*det)[2];\r
369 } else if(strcmp(opt, "pad")==0){\r
5591f3d6 370 Int_t stk(AliTRDgeometry::GetStack(idet));\r
371 xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
372 } else continue;\r
b9058d72 373 AliDebug(2, Form("det[%03d] 0[%+4.1f(%+4.1f) %+4.1f] 1[%+4.1f(%+4.1f) %+4.1f] opt[%d]", idet, xmin, (*det)[0], (*det)[1], xmax, (*det)[2], (*det)[3], iopt));\r
6687261d 374 if(iopt==1){\r
07b3a6a6 375 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
5591f3d6 376 gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
07b3a6a6 377 } else {\r
378 gdet->SetFillStyle(style[0]);\r
5591f3d6 379 gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
6687261d 380 if(iopt==2){\r
b9058d72 381 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
5591f3d6 382 gdet->DrawBox(xmin, (*det)[1], xmax, 0.5*((*det)[3]+(*det)[1]));\r
6687261d 383 } else if(iopt==3){\r
07b3a6a6 384 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
5591f3d6 385 gdet->DrawBox(xmin, 0.5*((*det)[3]+(*det)[1]), xmax, (*det)[3]);\r
07b3a6a6 386 } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
387 }\r
82e6e5dc 388 }\r
5591f3d6 389 Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
390 xmin=0.;\r
391 if(strcmp(opt, "pad")==0) xmin=38.;\r
392 TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
393 for(Int_t is(0); is<AliTRDgeometry::kNsector; is++) sm->DrawLatex(xmin, -TMath::Pi()+(is+0.5)*dsm, Form("%02d", is>=9?(is-9):(is+9)));\r
82e6e5dc 394}\r
395\r
396\r
397//_______________________________________________________\r
94b94be0 398void AliTRDrecoTask::MakeSummary()\r
399{\r
400// To be implemented by particular tasks\r
401 AliWarning("Summary not available");\r
402}\r
403\r
404//_______________________________________________________\r
405void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
406{\r
407// Generic debug handler\r
408\r
409 AliAnalysisTaskSE::SetDebugLevel(level);\r
410 if(DebugLevel()>=1){\r
411 AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
412 TDirectory *savedir = gDirectory;\r
413 fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
414 savedir->cd();\r
415 }\r
416}\r
417\r
418//____________________________________________________________________\r
419void AliTRDrecoTask::Terminate(Option_t *)\r
420{\r
421 //\r
422 // Terminate\r
423 //\r
424\r
425 if(fgDebugStream){ \r
426 delete fgDebugStream;\r
427 fgDebugStream = NULL;\r
428 }\r
429 fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
430 if(fContainer && fRunTerminate){\r
431 PostProcess();\r
432 MakeSummary();\r
433 }\r
434}\r
435\r
436//________________________________________________________\r
437void AliTRDrecoTask::Adjust(TF1 *f, TH1 * const h)\r
438{\r
439// Helper function to avoid duplication of code\r
440// Make first guesses on the fit parameters\r
441\r
442 // find the intial parameters of the fit !! (thanks George)\r
443 Int_t nbinsy = Int_t(.5*h->GetNbinsX());\r
444 Double_t sum = 0.;\r
445 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;\r
446 f->SetParLimits(0, 0., 3.*sum);\r
447 f->SetParameter(0, .9*sum);\r
448\r
449 f->SetParLimits(1, -.2, .2);\r
450 f->SetParameter(1, -0.1);\r
451\r
452 f->SetParLimits(2, 0., 4.e-1);\r
453 f->SetParameter(2, 2.e-2);\r
454 if(f->GetNpar() <= 4) return;\r
455\r
456 f->SetParLimits(3, 0., sum);\r
457 f->SetParameter(3, .1*sum);\r
458\r
459 f->SetParLimits(4, -.3, .3);\r
460 f->SetParameter(4, 0.);\r
461\r
462 f->SetParLimits(5, 0., 1.e2);\r
463 f->SetParameter(5, 2.e-1);\r
464}\r