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