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