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