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