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