]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TRD/AliTRDrecoTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Cleanup
[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
eb05d549 25#include "TH2.h"\r
26#include "TH3.h"\r
94b94be0 27#include "TF1.h"\r
28#include "TObjArray.h"\r
29#include "TDirectory.h"\r
30#include "TTreeStream.h"\r
82e6e5dc 31#include "TBox.h"\r
5591f3d6 32#include "TLatex.h"\r
82e6e5dc 33#include "TVectorT.h"\r
94b94be0 34\r
35#include "AliLog.h"\r
36#include "AliAnalysisTask.h"\r
100f4577 37#include "AliAnalysisManager.h"\r
3ceb45ae 38#include "AliExternalTrackParam.h"\r
94b94be0 39\r
40#include "info/AliTRDeventInfo.h"\r
847569f2 41#include "info/AliTRDtrendingManager.h"\r
94b94be0 42#include "AliTRDrecoTask.h"\r
3ceb45ae 43#include "AliTRDtrackV1.h"\r
44#include "AliTRDpidUtil.h"\r
94b94be0 45\r
46ClassImp(AliTRDrecoTask)\r
47\r
33056e04 48Float_t AliTRDrecoTask::fgPt0[AliTRDrecoTask::fgNPt0] = {0.5, 0.8, 1.5, 5};\r
94b94be0 49TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
50//_______________________________________________________\r
51AliTRDrecoTask::AliTRDrecoTask()\r
52 : AliAnalysisTaskSE()\r
53 ,fNRefFigures(0)\r
82e6e5dc 54 ,fDets(NULL)\r
33056e04 55 ,fDetsV(NULL)\r
94b94be0 56 ,fContainer(NULL)\r
57 ,fEvent(NULL)\r
58 ,fTracks(NULL)\r
f073d500 59 ,fClusters(NULL)\r
60 ,fkClusters(NULL)\r
94b94be0 61 ,fkTrack(NULL)\r
62 ,fkMC(NULL)\r
63 ,fkESD(NULL)\r
3ceb45ae 64 ,fSpecies(-6)\r
65 ,fPt(-1.)\r
66 ,fPhi(0.)\r
67 ,fEta(0.)\r
94b94be0 68 ,fPlotFuncList(NULL)\r
f073d500 69 ,fDetFuncList(NULL)\r
94b94be0 70 ,fRunTerminate(kFALSE)\r
71{\r
72// Default constructor\r
73 snprintf(fNameId, 10, "no name");\r
74}\r
75\r
76//_______________________________________________________\r
77AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
78 : AliAnalysisTaskSE(name)\r
79 ,fNRefFigures(0)\r
82e6e5dc 80 ,fDets(NULL)\r
33056e04 81 ,fDetsV(NULL)\r
94b94be0 82 ,fContainer(NULL)\r
83 ,fEvent(NULL)\r
84 ,fTracks(NULL)\r
f073d500 85 ,fClusters(NULL)\r
86 ,fkClusters(NULL)\r
94b94be0 87 ,fkTrack(NULL)\r
88 ,fkMC(NULL)\r
89 ,fkESD(NULL)\r
3ceb45ae 90 ,fSpecies(-6)\r
91 ,fPt(-1.)\r
92 ,fPhi(0.)\r
93 ,fEta(0.)\r
94b94be0 94 ,fPlotFuncList(NULL)\r
f073d500 95 ,fDetFuncList(NULL)\r
94b94be0 96 ,fRunTerminate(kFALSE)\r
97{\r
98// Constructor for all derived performance tasks\r
99\r
100 SetTitle(title);\r
101 snprintf(fNameId, 10, "no name");\r
102 DefineInput (1, TObjArray::Class()); // track list\r
103 DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
f073d500 104 DefineInput (3, TObjArray::Class()); // cluster list object\r
94b94be0 105 DefineOutput(1, TObjArray::Class()); // histogram list\r
106}\r
107\r
108//_______________________________________________________\r
109AliTRDrecoTask::~AliTRDrecoTask() \r
110{\r
111\r
112 // Generic task destructor\r
113\r
114 AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
115 if(fgDebugStream){ \r
116 delete fgDebugStream;\r
117 fgDebugStream = NULL;\r
118 }\r
119\r
120 if(fPlotFuncList){\r
121 fPlotFuncList->Delete();\r
122 delete fPlotFuncList;\r
123 fPlotFuncList = NULL;\r
124 }\r
f073d500 125 if(fDetFuncList){\r
126 fDetFuncList->Delete();\r
127 delete fDetFuncList;\r
128 fDetFuncList = NULL;\r
129 }\r
94b94be0 130 \r
82e6e5dc 131 if(fDets){\r
132 if(fDets->IsOwner()) fDets->Delete();\r
133 delete fDets;\r
134 fDets = NULL;\r
135 }\r
33056e04 136 if(fDetsV) delete fDetsV; fDetsV=NULL;\r
82e6e5dc 137\r
f0857a6a 138 if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
94b94be0 139 if(fContainer->IsOwner()) fContainer->Delete();\r
140 delete fContainer;\r
141 fContainer = NULL;\r
142 }\r
143\r
847569f2 144/* if(fgTrendPoint){\r
94b94be0 145 TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
146 fgTrendPoint->Write();\r
147 delete fgTrendPoint;\r
148 fgTrendPoint=NULL;\r
149 gFile->Close();\r
847569f2 150 }*/\r
94b94be0 151}\r
152\r
153//_______________________________________________________\r
154Int_t AliTRDrecoTask::GetNRefFigures() const \r
155{ \r
156 if(!fNRefFigures) AliWarning("No reference plots available.");\r
157 return fNRefFigures; \r
158} \r
159\r
33056e04 160//____________________________________________________________________\r
161Int_t AliTRDrecoTask::GetPtBinSignificant(Float_t pt)\r
162{\r
163// Get significant (very low, low, medium, high, very high) pt bin\r
164\r
165 Int_t ipt(0);\r
166 while(ipt<fgNPt0){\r
167 if(pt<fgPt0[ipt]) break;\r
168 ipt++;\r
169 }\r
170 return ipt-1;\r
171}\r
172\r
94b94be0 173//_______________________________________________________\r
174void AliTRDrecoTask::UserCreateOutputObjects()\r
175{\r
176 if(!HasFunctorList()) InitFunctorList();\r
177 fContainer = Histos();\r
178 PostData(1, fContainer);\r
179}\r
180\r
181//_______________________________________________________\r
182void AliTRDrecoTask::UserExec(Option_t *)\r
183{\r
184// Loop over Plot functors published by particular tasks\r
185\r
f073d500 186 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));\r
187 fEvent = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
188 fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
94b94be0 189\r
190 if(!fPlotFuncList){\r
f073d500 191 AliWarning("No track functor list defined for the task");\r
94b94be0 192 return;\r
193 }\r
194 if(!fTracks) return;\r
195 if(!fTracks->GetEntriesFast()) return;\r
196 else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
197\r
178d284a 198 Int_t itrk(-1);\r
199 AliTRDtrackInfo *trackInfo(NULL);\r
94b94be0 200 TIter plotIter(fPlotFuncList);\r
201 TObjArrayIter trackIter(fTracks);\r
202 while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
178d284a 203 itrk++; fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
94b94be0 204 fkMC = trackInfo->GetMCinfo();\r
205 fkESD = trackInfo->GetESDinfo();\r
178d284a 206 if((fkTrack = trackInfo->GetTrack())){\r
207 // cache properties of the track at TRD entrance\r
208 // check input track status\r
209 AliExternalTrackParam *tin(NULL);\r
210 if(!(tin = fkTrack->GetTrackIn())) AliDebug(2, Form("Missing TRD track[%d] :: entry point.", itrk));\r
211 else {\r
212 fPt = tin->Pt();\r
213 fEta = tin->Eta();\r
214 Double_t xyz[3];\r
215 if(!tin->GetXYZ(xyz)) AliDebug(2, Form("Failed TRD track[%d] :: global track postion", itrk));\r
216 else fPhi = TMath::ATan2(xyz[1], xyz[0]);\r
217 fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
218 }\r
219 } else AliDebug(2, Form("Missing TRD track[%d].", itrk));\r
220\r
221 TMethodCall *plot(NULL);\r
94b94be0 222 plotIter.Reset();\r
178d284a 223 while((plot=dynamic_cast<TMethodCall*>(plotIter()))) plot->Execute(this);\r
94b94be0 224 }\r
f073d500 225 if(!fClusters) return;\r
226 if(!fDetFuncList){\r
178d284a 227 AliDebug(1, "No detector functor list defined for task");\r
f073d500 228 return;\r
229 }\r
230 TIter detIter(fDetFuncList);\r
231 for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
232 if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
233 TMethodCall *det(NULL);\r
234 detIter.Reset();\r
178d284a 235 while((det=dynamic_cast<TMethodCall*>(detIter()))) det->Execute(this);\r
f073d500 236 }\r
94b94be0 237}\r
238\r
239//_______________________________________________________\r
240Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
241{\r
242 AliWarning("Retrieving reference figures not implemented.");\r
243 return kFALSE;\r
244}\r
245\r
246//_______________________________________________________\r
847569f2 247Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val, Double_t err)\r
94b94be0 248{\r
249// Generic publisher for trend values\r
250\r
847569f2 251 AliTRDtrendingManager *tm = AliTRDtrendingManager::Instance();\r
252 tm->AddValue(Form("%s_%s", GetName(), name), val, err);\r
94b94be0 253 return kTRUE;\r
254}\r
255\r
256//_______________________________________________________\r
257void AliTRDrecoTask::InitFunctorList()\r
258{\r
259// Initialize list of functors\r
260\r
261 TClass *c = this->IsA();\r
3ed01fbe 262 if(fPlotFuncList) fPlotFuncList->Clear();\r
f073d500 263 if(fDetFuncList) fDetFuncList->Clear();\r
94b94be0 264\r
f073d500 265 TMethod *m(NULL);\r
94b94be0 266 TIter methIter(c->GetListOfMethods());\r
267 while((m=dynamic_cast<TMethod*>(methIter()))){\r
268 TString name(m->GetName());\r
f073d500 269 if(name.BeginsWith("Plot")){\r
270 if(!fPlotFuncList) fPlotFuncList = new TList();\r
271 fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
272 } else if(name.BeginsWith("Det")){\r
273 if(!fDetFuncList) fDetFuncList = new TList();\r
274 fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
275 }\r
94b94be0 276 }\r
277}\r
278\r
279//_______________________________________________________\r
280Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
281{\r
282// Generic container loader\r
283\r
284 if(!TFile::Open(file)){\r
285 AliWarning(Form("Couldn't open file %s.", file));\r
286 return kFALSE;\r
287 }\r
288 if(!gFile->cd(dir)){\r
289 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
290 return kFALSE;\r
291 }\r
292 TObjArray *o = NULL;\r
293 if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
294 AliWarning("Missing histogram container.");\r
295 return kFALSE;\r
296 }\r
297 fContainer = (TObjArray*)o->Clone(GetName());\r
298 gFile->Close();\r
299 return kTRUE;\r
300}\r
301\r
82e6e5dc 302//________________________________________________________\r
303Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
304{\r
305// Load detector map.\r
306\r
307 if(!TFile::Open(file)){\r
308 AliWarning(Form("Couldn't open file %s.", file));\r
309 return kFALSE;\r
310 }\r
311 if(!gFile->cd(dir)){\r
312 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
313 return kFALSE;\r
314 }\r
315 TObjArray *info = NULL;\r
316 if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
317 AliWarning("Missing TRDinfoGen container.");\r
318 return kFALSE;\r
319 }\r
320 TObjArray *dets = (TObjArray*)info->FindObject("Chambers");\r
321 if(!dets){\r
33056e04 322 TVector *vdets = (TVector*)info->At(4);\r
323 if(!vdets){\r
324 AliWarning("Missing detector map from TRDinfoGen results.");\r
325 return kFALSE;\r
326 } else fDetsV = (TVector*)vdets->Clone();\r
327 } else fDets = (TObjArray*)dets->Clone("Chambers");\r
82e6e5dc 328 gFile->Close();\r
329 return kTRUE;\r
330}\r
331\r
332\r
94b94be0 333//________________________________________________________\r
334Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
335 //\r
336 // Store the output graphs in a ROOT file\r
337 // Input TObject array will not be written as Key to the file,\r
338 // only content itself\r
339 //\r
340\r
341 TDirectory *cwd = gDirectory;\r
342 if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
343\r
344 TIterator *iter = results->MakeIterator();\r
345 TObject *inObject = NULL, *outObject = NULL;\r
346 while((inObject = iter->Next())){\r
347 outObject = inObject->Clone();\r
348 outObject->Write(NULL, TObject::kSingleKey);\r
349 }\r
350 delete iter;\r
351 gFile->Close(); delete gFile;\r
352 cwd->cd(); \r
353 return kTRUE;\r
354}\r
355\r
356//_______________________________________________________\r
357Bool_t AliTRDrecoTask::PostProcess()\r
358{\r
359// To be implemented by particular tasks\r
360\r
361 AliWarning("Post processing of reference histograms not implemented.");\r
362 return kTRUE;\r
363}\r
364\r
82e6e5dc 365//_______________________________________________________\r
5591f3d6 366void AliTRDrecoTask::MakeDetectorPlot(Int_t ly, const Option_t *opt)\r
82e6e5dc 367{\r
368// Draw chamber boundaries in eta/phi plots with misalignments\r
369// based on info collected by AliTRDinfoGen\r
370\r
371 if(!fDets){\r
33056e04 372 AliWarning("OLD Detector map and status not available. Try NEW");\r
373 MakeDetectorPlotNEW(ly, opt);\r
82e6e5dc 374 return;\r
375 }\r
5591f3d6 376 Float_t xmin(0.), xmax(0.);\r
82e6e5dc 377 TBox *gdet = new TBox();\r
eb05d549 378 gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
82e6e5dc 379 Int_t style[] = {0, 3003};\r
380 for(Int_t idet(0); idet<540; idet++){\r
381 if(idet%6 != ly) continue;\r
382 TVectorF *det((TVectorF*)fDets->At(idet));\r
383 if(!det) continue;\r
82e6e5dc 384 Int_t iopt = Int_t((*det)[4]);\r
b9058d72 385 if(strcmp(opt, "eta")==0){\r
386 xmin=(*det)[0]; xmax=(*det)[2];\r
387 } else if(strcmp(opt, "pad")==0){\r
5591f3d6 388 Int_t stk(AliTRDgeometry::GetStack(idet));\r
389 xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
390 } else continue;\r
b9058d72 391 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 392 if(iopt==1){\r
07b3a6a6 393 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
5591f3d6 394 gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
07b3a6a6 395 } else {\r
396 gdet->SetFillStyle(style[0]);\r
5591f3d6 397 gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
6687261d 398 if(iopt==2){\r
b9058d72 399 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
5591f3d6 400 gdet->DrawBox(xmin, (*det)[1], xmax, 0.5*((*det)[3]+(*det)[1]));\r
6687261d 401 } else if(iopt==3){\r
07b3a6a6 402 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
5591f3d6 403 gdet->DrawBox(xmin, 0.5*((*det)[3]+(*det)[1]), xmax, (*det)[3]);\r
07b3a6a6 404 } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
405 }\r
82e6e5dc 406 }\r
5591f3d6 407 Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
408 xmin=0.;\r
409 if(strcmp(opt, "pad")==0) xmin=38.;\r
410 TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
411 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 412}\r
413\r
33056e04 414//_______________________________________________________\r
415void AliTRDrecoTask::MakeDetectorPlotNEW(Int_t ly, const Option_t *opt)\r
416{\r
417// Draw chamber boundaries in eta/phi plots with misalignments\r
418// based on info collected by AliTRDinfoGen NEW data storage\r
419\r
420 if(!fDetsV){\r
421 AliWarning("NEW Detector map and status not available.");\r
422 return;\r
423 }\r
424 Float_t xmin(0.), xmax(0.);\r
425 TBox *gdet = new TBox();\r
426 gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
427 Int_t style[] = {0, 3003};\r
428 for(Int_t idet(0), jdet(0); idet<AliTRDgeometry::kNdet; idet++, jdet+=5){\r
429 if(idet%6 != ly) continue;\r
430 Int_t iopt = Int_t((*fDetsV)[jdet+4]);\r
431 if(strcmp(opt, "eta")==0){\r
432 xmin=(*fDetsV)[jdet+0]; xmax=(*fDetsV)[jdet+2];\r
433 } else if(strcmp(opt, "pad")==0){\r
434 Int_t stk(AliTRDgeometry::GetStack(idet));\r
435 xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
436 } else continue;\r
437 AliDebug(2, Form("det[%03d] 0[%+4.1f(%+4.1f) %+4.1f] 1[%+4.1f(%+4.1f) %+4.1f] opt[%d]", idet, xmin, (*fDetsV)[jdet+0], (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+2], (*fDetsV)[jdet+3], iopt));\r
438 if(iopt==1){\r
439 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
440 gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+3]);\r
441 } else {\r
442 gdet->SetFillStyle(style[0]);\r
443 gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+3]);\r
444 if(iopt==2){\r
445 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
446 gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, 0.5*((*fDetsV)[jdet+3]+(*fDetsV)[jdet+1]));\r
447 } else if(iopt==3){\r
448 gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
449 gdet->DrawBox(xmin, 0.5*((*fDetsV)[jdet+3]+(*fDetsV)[jdet+1]), xmax, (*fDetsV)[jdet+3]);\r
450 } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
451 }\r
452 }\r
453 Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
454 xmin=0.;\r
455 if(strcmp(opt, "pad")==0) xmin=38.;\r
456 TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
457 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
458}\r
459\r
82e6e5dc 460\r
94b94be0 461//_______________________________________________________\r
462void AliTRDrecoTask::MakeSummary()\r
463{\r
464// To be implemented by particular tasks\r
465 AliWarning("Summary not available");\r
466}\r
467\r
468//_______________________________________________________\r
469void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
470{\r
471// Generic debug handler\r
472\r
473 AliAnalysisTaskSE::SetDebugLevel(level);\r
474 if(DebugLevel()>=1){\r
475 AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
476 TDirectory *savedir = gDirectory;\r
477 fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
478 savedir->cd();\r
479 }\r
480}\r
481\r
482//____________________________________________________________________\r
483void AliTRDrecoTask::Terminate(Option_t *)\r
484{\r
485 //\r
486 // Terminate\r
487 //\r
488\r
489 if(fgDebugStream){ \r
490 delete fgDebugStream;\r
491 fgDebugStream = NULL;\r
492 }\r
493 fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
494 if(fContainer && fRunTerminate){\r
495 PostProcess();\r
496 MakeSummary();\r
497 }\r
498}\r
499\r
eb05d549 500//________________________________________________________\r
cc98ff07 501Float_t AliTRDrecoTask::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)\r
eb05d549 502{\r
503// Normalize histo content to the mean value in the range specified by bin ranges\r
504// [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.\r
505// Optionally a threshold "thr" can be specified to disregard entries with no meaning\r
506\r
507 Float_t s = 0., c=0.; Int_t is(0);\r
508 for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){\r
509 for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){\r
510 if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
511 s += c; is++;\r
512 }\r
513 }\r
ea90ffcd 514 s/= (is?is:1);\r
eb05d549 515 for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
516 for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
ea90ffcd 517 if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-1000);\r
eb05d549 518 else h2->SetBinContent(ix, iy, 100.*(c/s-1.));\r
519 }\r
520 }\r
cc98ff07 521 return s;\r
eb05d549 522}\r
523\r
524//________________________________________________________\r
525void AliTRDrecoTask::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)\r
526{\r
527// Set range on Z axis such to avoid outliers\r
528\r
529 Float_t c(0.), dz(1.e-3*(max-min));\r
530 for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
531 for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
532 if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
533 if(c<=min) h2->SetBinContent(ix, iy, min+dz);\r
534 }\r
535 }\r
536 h2->GetZaxis()->SetRangeUser(min, max);\r
537}\r
538\r
539//________________________________________________________\r
540Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)\r
541{\r
542// return mean number of entries/bin of histogram "h"\r
543// if option "opt" is given the following values are accepted:\r
544// "<" : consider only entries less than "cut"\r
545// ">" : consider only entries greater than "cut"\r
546\r
547 //Int_t dim(h->GetDimension());\r
548 Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());\r
549 Double_t sum(0.); Int_t n(0);\r
550 for(Int_t ix(1); ix<=nbx; ix++)\r
551 for(Int_t iy(1); iy<=nby; iy++)\r
552 for(Int_t iz(1); iz<=nbz; iz++){\r
553 if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}\r
554 else{\r
555 if(strcmp(opt, "<")==0) {\r
556 if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
557 } else if(strcmp(opt, ">")==0){\r
558 if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
559 } else {sum += h->GetBinContent(ix, iy, iz); n++;}\r
560 }\r
561 }\r
562 return n>0?sum/n:0.;\r
563}\r
564\r
565//________________________________________________________\r
566AliTRDrecoTask::AliTRDrecoProjection::AliTRDrecoProjection()\r
567 :TNamed()\r
568 ,fH(NULL)\r
569 ,fNrebin(0)\r
570 ,fRebinX(NULL)\r
571 ,fRebinY(NULL)\r
572{\r
573 // constructor\r
574 memset(fAx, 0, 3*sizeof(Int_t));\r
575 memset(fRange, 0, 4*sizeof(Float_t));\r
576}\r
577\r
578//________________________________________________________\r
579AliTRDrecoTask::AliTRDrecoProjection::~AliTRDrecoProjection()\r
580{\r
581 // destructor\r
582 if(fH) delete fH;\r
583}\r
584\r
585//________________________________________________________\r
586void AliTRDrecoTask::AliTRDrecoProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])\r
587{\r
588// check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"\r
589 if(!aa[ix] || !aa[iy] || !aa[iz]) return;\r
590 TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);\r
591 // check ax definiton to protect against older versions of the data\r
592 if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){\r
593 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));\r
594 return;\r
595 }\r
596 if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){\r
597 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));\r
598 return;\r
599 }\r
600 if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){\r
601 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
602 return;\r
603 }\r
604 SetNameTitle(n,t);\r
605 fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),\r
606 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
607 ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
608 az->GetNbins(), az->GetXmin(), az->GetXmax());\r
609 fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;\r
610 fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;\r
611 AliDebug(2, Form("H3(%s, %s) :: %s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]", n, t,\r
612 ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
613 ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
614 az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
615}\r
616\r
617//________________________________________________________\r
618AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator=(const AliTRDrecoProjection& rhs)\r
619{\r
620// copy projections\r
621 if(this == &rhs) return *this;\r
622\r
623 TNamed::operator=(rhs);\r
624 if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}\r
625 if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);\r
626 memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));\r
627 memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));\r
628 if(fH) delete fH;\r
629 if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));\r
630 return *this;\r
631}\r
632\r
633//________________________________________________________\r
634AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator+=(const AliTRDrecoProjection& other)\r
635{\r
636// increment projections\r
637 if(!fH || !other.fH) return *this;\r
638 AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));\r
639 fH->Add(other.fH);\r
640 return *this;\r
641}\r
642\r
643//________________________________________________________\r
644void AliTRDrecoTask::AliTRDrecoProjection::Increment(Int_t bin[], Double_t v)\r
645{\r
646// increment bin with value "v" pointed by general coord in "bin"\r
647 if(!fH) return;\r
648 AliDebug(4, Form(" %s[%2d]", fH->GetName(), Int_t(v)));\r
649 fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));\r
650}\r
651\r
668a0654 652//________________________________________________________\r
653Double_t AliTRDrecoTask::AliTRDrecoProjection::GetTrendValue(const Int_t mid, Double_t *m, Double_t *s) const\r
654{\r
655// Return result of fitting the main distribution (represented on the z axis) with the function selected\r
656// "mid". Optionally return the Mean and RMS of the distribution pointing to "m" and "s"\r
657\r
658 if(!fH){\r
659 AliDebug(1, Form("Missing 3D in %s", GetName()));\r
660 return -999.;\r
661 }\r
662 TH1 *h1s(NULL);\r
f036f7ec 663 if(!(h1s = (TH1D*)fH->Project3D("z"))){\r
668a0654 664 AliDebug(1, Form("Failed Project3D(\"z\") in %s", GetName()));\r
665 return -999.;\r
666 }\r
667 Int_t ne((Int_t)h1s->Integral());\r
f036f7ec 668 if(ne<30){\r
668a0654 669 AliDebug(1, Form("Statistics too low[%2d] in %s", ne, GetName()));\r
670 return -999.;\r
671 }\r
672 TAxis *az(h1s->GetXaxis());\r
673 Float_t vm(h1s->GetMean()), v(vm), ve(h1s->GetRMS());\r
674 if(mid==1){\r
675 TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
676 fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, vm); fg.SetParameter(2, ve);\r
677 h1s->Fit(&fg, "WQ0");\r
678 v = fg.GetParameter(1);\r
679 } else if (mid==2) {\r
680 TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
681 fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, vm); fl.SetParameter(2, ve);\r
682 h1s->Fit(&fl, "WQ0");\r
683 v = fl.GetMaximumX();\r
684 }\r
685 if(m) *m = vm;\r
686 if(s) *s = ve;\r
f036f7ec 687 AliDebug(2, Form("%s[%d]:: %f {%f %f} Entries[%d]", fH->GetName(), mid, v, m?(*m):0., s?(*s):0., (Int_t)h1s->Integral()));\r
688\r
668a0654 689 return v;\r
690}\r
691\r
eb05d549 692//________________________________________________________\r
693TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)\r
694{\r
695// build the 2D projection and adjust binning\r
696\r
697 const Char_t *title[] = {"Mean", "#mu", "MPV"};\r
ea90ffcd 698 if(!fH){\r
699 AliDebug(1, Form("Missing 3D in %s", GetName()));\r
700 return NULL;\r
701 }\r
eb05d549 702 TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());\r
703 TH2D *h2s(NULL), *hyx(NULL);\r
ea90ffcd 704 if(!(h2s = (TH2D*)fH->Project3D("yx"))){\r
705 AliDebug(1, Form("Failed Project3D(\"yx\") in %s", GetName()));\r
706 return NULL;\r
707 }\r
eb05d549 708 // save a copy of the original distribution\r
709 if(!del){\r
710 hyx = (TH2D*)h2s->Clone();\r
711 hyx->SetName(Form("%sEn", fH->GetName()));\r
712 }\r
713 Int_t irebin(0), dxBin(1), dyBin(1);\r
714 while(irebin<fNrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, ">")<nstat)){\r
715 h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);\r
716 dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];\r
717 irebin++;\r
718 }\r
719 Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());\r
720 delete h2s;\r
721 if(mid<0) return NULL;\r
722\r
723 // start projection\r
724 TH1 *h(NULL); Int_t n(0);\r
725 Float_t dz=(fRange[1]-fRange[1])/ncol;\r
726 TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));\r
727 Int_t nt(tokenTitle->GetEntriesFast());\r
f985021d 728 TH2 *h2(NULL);\r
729 if((h2 = (TH2*)gDirectory->Get(Form("%s_2D", fH->GetName())))) delete h2; // avoid ROOT warning messages\r
730 h2 = new TH2F(Form("%s_2D", fH->GetName()),\r
eb05d549 731 Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], nt>0?(*tokenTitle)[0]->GetName():"", nt>1?(*tokenTitle)[1]->GetName():""),\r
732 nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());\r
733 h2->SetContour(ncol);\r
734 h2->GetZaxis()->CenterTitle();\r
735 h2->GetZaxis()->SetTitleOffset(1.4);\r
736 h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);\r
737 AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));\r
738 for(Int_t iy(0); iy<ny; iy++){\r
739 for(Int_t ix(0); ix<nx; ix++){\r
ea90ffcd 740 h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin);\r
eb05d549 741 Int_t ne((Int_t)h->Integral());\r
ea90ffcd 742 //printf(" x[%2d %2d] y[%2d %2d] ne[%4d]\n", ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin, ne);\r
eb05d549 743 if(ne<nstat/2){\r
744 h2->SetBinContent(ix+1, iy+1, -999);\r
745 h2->SetBinError(ix+1, iy+1, 1.);\r
746 n++;\r
747 }else{\r
ea90ffcd 748 // redo the projection by adding 1 bin @ left and 1 bin @ right for smoothing\r
749 h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin, (ix+1)*dxBin+1, iy*dyBin, (iy+1)*dyBin+1);\r
eb05d549 750 Float_t v(h->GetMean()), ve(h->GetRMS());\r
751 if(mid==1){\r
752 TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
753 fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);\r
33056e04 754 h->Fit(&fg, "WQ0");\r
eb05d549 755 v = fg.GetParameter(1); ve = fg.GetParameter(2);\r
756 } else if (mid==2) {\r
757 TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
758 fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);\r
33056e04 759 h->Fit(&fl, "WQ0");\r
eb05d549 760 v = fl.GetMaximumX(); ve = fl.GetParameter(2);\r
761/* TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());\r
762 fgle.SetParameter(0, fl.GetParameter(0));\r
763 fgle.SetParameter(1, fl.GetParameter(1));\r
764 fgle.SetParameter(2, fl.GetParameter(2));\r
765 fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);\r
766 h->Fit(&fgle, "WQ");\r
767 v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/\r
768 }\r
769 if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);\r
770 else h2->SetBinContent(ix+1, iy+1, v);\r
771 h2->SetBinError(ix+1, iy+1, ve);\r
772 }\r
773 }\r
774 }\r
775 if(h) delete h;\r
ea90ffcd 776 if(n==nx*ny){ // clean empty projections\r
777 AliDebug(1, Form("Empty projection in %s", GetName()));\r
778 delete h2; h2=NULL;\r
779 }\r
eb05d549 780 return h2;\r
781}\r
782\r
783//________________________________________________________\r
784void AliTRDrecoTask::AliTRDrecoProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])\r
785{\r
786// define rebinning strategy for this projection\r
787 fNrebin = n;\r
788 fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));\r
789 fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));\r
790}\r
791\r
792\r
793\r