add cluster monitoring from RecPoints
[u/mrichter/AliRoot.git] / PWG1 / 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 "TF1.h"\r
26 #include "TObjArray.h"\r
27 #include "TDirectory.h"\r
28 #include "TTreeStream.h"\r
29 #include "TBox.h"\r
30 #include "TVectorT.h"\r
31 \r
32 #include "AliLog.h"\r
33 #include "AliAnalysisTask.h"\r
34 #include "AliAnalysisManager.h"\r
35 #include "AliExternalTrackParam.h"\r
36 \r
37 #include "info/AliTRDeventInfo.h"\r
38 #include "AliTRDrecoTask.h"\r
39 #include "AliTRDtrackV1.h"\r
40 #include "AliTRDpidUtil.h"\r
41 \r
42 ClassImp(AliTRDrecoTask)\r
43 \r
44 TList* AliTRDrecoTask::fgTrendPoint(NULL);\r
45 TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
46 //_______________________________________________________\r
47 AliTRDrecoTask::AliTRDrecoTask()\r
48   : AliAnalysisTaskSE()\r
49   ,fNRefFigures(0)\r
50   ,fDets(NULL)\r
51   ,fContainer(NULL)\r
52   ,fEvent(NULL)\r
53   ,fTracks(NULL)\r
54   ,fClusters(NULL)\r
55   ,fkClusters(NULL)\r
56   ,fkTrack(NULL)\r
57   ,fkMC(NULL)\r
58   ,fkESD(NULL)\r
59   ,fSpecies(-6)\r
60   ,fPt(-1.)\r
61   ,fPhi(0.)\r
62   ,fEta(0.)\r
63   ,fPlotFuncList(NULL)\r
64   ,fDetFuncList(NULL)\r
65   ,fRunTerminate(kFALSE)\r
66 {\r
67 // Default constructor\r
68   snprintf(fNameId, 10, "no name");\r
69 }\r
70 \r
71 //_______________________________________________________\r
72 AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
73   : AliAnalysisTaskSE(name)\r
74   ,fNRefFigures(0)\r
75   ,fDets(NULL)\r
76   ,fContainer(NULL)\r
77   ,fEvent(NULL)\r
78   ,fTracks(NULL)\r
79   ,fClusters(NULL)\r
80   ,fkClusters(NULL)\r
81   ,fkTrack(NULL)\r
82   ,fkMC(NULL)\r
83   ,fkESD(NULL)\r
84   ,fSpecies(-6)\r
85   ,fPt(-1.)\r
86   ,fPhi(0.)\r
87   ,fEta(0.)\r
88   ,fPlotFuncList(NULL)\r
89   ,fDetFuncList(NULL)\r
90   ,fRunTerminate(kFALSE)\r
91 {\r
92 // Constructor for all derived performance tasks\r
93 \r
94   SetTitle(title);\r
95   snprintf(fNameId, 10, "no name");\r
96   DefineInput (1, TObjArray::Class()); // track list\r
97   DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
98   DefineInput (3, TObjArray::Class()); // cluster list object\r
99   DefineOutput(1, TObjArray::Class()); // histogram list\r
100 }\r
101 \r
102 //_______________________________________________________\r
103 AliTRDrecoTask::~AliTRDrecoTask() \r
104 {\r
105 \r
106   // Generic task destructor\r
107 \r
108   AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
109   if(fgDebugStream){ \r
110     delete fgDebugStream;\r
111     fgDebugStream = NULL;\r
112   }\r
113 \r
114   if(fPlotFuncList){\r
115     fPlotFuncList->Delete();\r
116     delete fPlotFuncList;\r
117     fPlotFuncList = NULL;\r
118   }\r
119   if(fDetFuncList){\r
120     fDetFuncList->Delete();\r
121     delete fDetFuncList;\r
122     fDetFuncList = NULL;\r
123   }\r
124   \r
125   if(fDets){\r
126     if(fDets->IsOwner()) fDets->Delete();\r
127     delete fDets;\r
128     fDets = NULL;\r
129   }\r
130 \r
131   if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
132     if(fContainer->IsOwner()) fContainer->Delete();\r
133     delete fContainer;\r
134     fContainer = NULL;\r
135   }\r
136 \r
137   if(fgTrendPoint){\r
138     TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
139     fgTrendPoint->Write();\r
140     delete fgTrendPoint;\r
141     fgTrendPoint=NULL;\r
142     gFile->Close();\r
143   }\r
144 }\r
145 \r
146 //_______________________________________________________\r
147 Int_t AliTRDrecoTask::GetNRefFigures() const  \r
148\r
149   if(!fNRefFigures) AliWarning("No reference plots available.");\r
150   return fNRefFigures; \r
151\r
152 \r
153 //_______________________________________________________\r
154 void AliTRDrecoTask::UserCreateOutputObjects()\r
155 {\r
156   if(!HasFunctorList()) InitFunctorList();\r
157   fContainer = Histos();\r
158   PostData(1, fContainer);\r
159 }\r
160 \r
161 //_______________________________________________________\r
162 void AliTRDrecoTask::UserExec(Option_t *)\r
163 {\r
164 // Loop over Plot functors published by particular tasks\r
165 \r
166   fTracks   = dynamic_cast<TObjArray *>(GetInputData(1));\r
167   fEvent    = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
168   fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
169 \r
170   if(!fPlotFuncList){\r
171     AliWarning("No track functor list defined for the task");\r
172     return;\r
173   }\r
174   if(!fTracks) return;\r
175   if(!fTracks->GetEntriesFast()) return;\r
176   else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
177 \r
178   AliTRDtrackInfo *trackInfo = NULL;\r
179   TIter plotIter(fPlotFuncList);\r
180   TObjArrayIter trackIter(fTracks);\r
181   while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
182     fkTrack = trackInfo->GetTrack();\r
183     fkMC    = trackInfo->GetMCinfo();\r
184     fkESD   = trackInfo->GetESDinfo();\r
185     // cache properties of the track at TRD entrance\r
186     // check input track status\r
187     AliExternalTrackParam *tin(NULL);\r
188     fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
189     if(!fkTrack || !(tin = fkTrack->GetTrackIn())) AliDebug(2, "Track did not entered TRD fiducial volume.");\r
190     else {\r
191       fPt   = tin->Pt();\r
192       fEta  = tin->Eta();\r
193       Double_t xyz[3];\r
194       if(!tin->GetXYZ(xyz)) AliDebug(2, "Failed getting global track postion");\r
195       else fPhi  = TMath::ATan2(xyz[1], xyz[0]);\r
196       fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
197     }\r
198 \r
199     TMethodCall *plot = NULL;\r
200     plotIter.Reset();\r
201     while((plot=dynamic_cast<TMethodCall*>(plotIter()))){\r
202       plot->Execute(this);\r
203     }\r
204   }\r
205   if(!fClusters) return;\r
206   if(!fDetFuncList){\r
207     AliDebug(1, "No detector functor list defined for the task");\r
208     return;\r
209   }\r
210   TIter detIter(fDetFuncList);\r
211   for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
212     if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
213     TMethodCall *det(NULL);\r
214     detIter.Reset();\r
215     while((det=dynamic_cast<TMethodCall*>(detIter()))){\r
216       det->Execute(this);\r
217     }\r
218   }\r
219 }\r
220 \r
221 //_______________________________________________________\r
222 Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
223 {\r
224   AliWarning("Retrieving reference figures not implemented.");\r
225   return kFALSE;\r
226 }\r
227 \r
228 //_______________________________________________________\r
229 Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val)\r
230 {\r
231 // Generic publisher for trend values\r
232 \r
233   if(!fgTrendPoint){\r
234     fgTrendPoint = new TList();\r
235     fgTrendPoint->SetOwner();\r
236   }\r
237   fgTrendPoint->AddLast(new TNamed(Form("%s_%s", GetName(), name), Form("%f", val)));\r
238   return kTRUE;\r
239 }\r
240 \r
241 //_______________________________________________________\r
242 void AliTRDrecoTask::InitFunctorList()\r
243 {\r
244 // Initialize list of functors\r
245 \r
246   TClass *c = this->IsA();\r
247   if(fPlotFuncList) fPlotFuncList->Clear();\r
248   if(fDetFuncList) fDetFuncList->Clear();\r
249 \r
250   TMethod *m(NULL);\r
251   TIter methIter(c->GetListOfMethods());\r
252   while((m=dynamic_cast<TMethod*>(methIter()))){\r
253     TString name(m->GetName());\r
254     if(name.BeginsWith("Plot")){\r
255       if(!fPlotFuncList) fPlotFuncList = new TList();\r
256       fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
257     } else if(name.BeginsWith("Det")){\r
258       if(!fDetFuncList) fDetFuncList = new TList();\r
259       fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
260     }\r
261   }\r
262 }\r
263 \r
264 //_______________________________________________________\r
265 Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
266 {\r
267 // Generic container loader\r
268 \r
269   if(!TFile::Open(file)){\r
270     AliWarning(Form("Couldn't open file %s.", file));\r
271     return kFALSE;\r
272   }\r
273   if(!gFile->cd(dir)){\r
274     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
275     return kFALSE;\r
276   }\r
277   TObjArray *o = NULL;\r
278   if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
279     AliWarning("Missing histogram container.");\r
280     return kFALSE;\r
281   }\r
282   fContainer = (TObjArray*)o->Clone(GetName());\r
283   gFile->Close();\r
284   return kTRUE;\r
285 }\r
286 \r
287 //________________________________________________________\r
288 Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
289 {\r
290 // Load detector map.\r
291 \r
292   if(!TFile::Open(file)){\r
293     AliWarning(Form("Couldn't open file %s.", file));\r
294     return kFALSE;\r
295   }\r
296   if(!gFile->cd(dir)){\r
297     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
298     return kFALSE;\r
299   }\r
300   TObjArray *info = NULL;\r
301   if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
302     AliWarning("Missing TRDinfoGen container.");\r
303     return kFALSE;\r
304   }\r
305   TObjArray *dets = (TObjArray*)info->FindObject("Chambers");\r
306   if(!dets){\r
307     AliWarning("Missing detector map from TRDinfoGen results.");\r
308     info->ls();\r
309     return kFALSE;\r
310   }\r
311   fDets = (TObjArray*)dets->Clone("Chambers");\r
312   gFile->Close();\r
313   return kTRUE;\r
314 }\r
315 \r
316 \r
317 //________________________________________________________\r
318 Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
319   //\r
320   // Store the output graphs in a ROOT file\r
321   // Input TObject array will not be written as Key to the file,\r
322   // only content itself\r
323   //\r
324 \r
325   TDirectory *cwd = gDirectory;\r
326   if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
327 \r
328   TIterator *iter = results->MakeIterator();\r
329   TObject *inObject = NULL, *outObject = NULL;\r
330   while((inObject = iter->Next())){\r
331     outObject = inObject->Clone();\r
332     outObject->Write(NULL, TObject::kSingleKey);\r
333   }\r
334   delete iter;\r
335   gFile->Close(); delete gFile;\r
336   cwd->cd(); \r
337   return kTRUE;\r
338 }\r
339 \r
340 //_______________________________________________________\r
341 Bool_t AliTRDrecoTask::PostProcess()\r
342 {\r
343 // To be implemented by particular tasks\r
344 \r
345   AliWarning("Post processing of reference histograms not implemented.");\r
346   return kTRUE;\r
347 }\r
348 \r
349 //_______________________________________________________\r
350 void AliTRDrecoTask::MakeDetectorPlot(Int_t ly)\r
351 {\r
352 // Draw chamber boundaries in eta/phi plots with misalignments\r
353 // based on info collected by AliTRDinfoGen\r
354 \r
355   if(!fDets){\r
356     AliWarning("Detector map and status not available.");\r
357     return;\r
358   }\r
359 \r
360   TBox *gdet = new TBox();\r
361   gdet->SetLineWidth(kBlack);gdet->SetFillColor(kBlack);\r
362   Int_t style[] = {0, 3003};\r
363   for(Int_t idet(0); idet<540; idet++){\r
364     if(idet%6 != ly) continue;\r
365     TVectorF *det((TVectorF*)fDets->At(idet));\r
366     if(!det) continue;\r
367     Int_t iopt = Int_t((*det)[4]);\r
368     AliDebug(2, Form("det[%03d] 0[%+4.1f %+4.1f] 1[%+4.1f %+4.1f] opt[%d]", idet, (*det)[0], (*det)[1], (*det)[2], (*det)[3], iopt));\r
369     if(iopt==1){\r
370       gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
371       gdet->DrawBox((*det)[0], (*det)[1], (*det)[2], (*det)[3]);\r
372     } else {\r
373       gdet->SetFillStyle(style[0]);\r
374       gdet->DrawBox((*det)[0], (*det)[1], (*det)[2], (*det)[3]);\r
375       if(iopt==2){\r
376         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlue);\r
377         gdet->DrawBox((*det)[0], (*det)[1], (*det)[2], 0.5*((*det)[3]+(*det)[1]));\r
378       } else if(iopt==3){\r
379         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
380         gdet->DrawBox((*det)[0], 0.5*((*det)[3]+(*det)[1]), (*det)[2], (*det)[3]);\r
381       } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
382     }\r
383   }\r
384 }\r
385 \r
386 \r
387 //_______________________________________________________\r
388 void AliTRDrecoTask::MakeSummary()\r
389 {\r
390 // To be implemented by particular tasks\r
391   AliWarning("Summary not available");\r
392 }\r
393 \r
394 //_______________________________________________________\r
395 void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
396 {\r
397 // Generic debug handler\r
398 \r
399   AliAnalysisTaskSE::SetDebugLevel(level);\r
400   if(DebugLevel()>=1){\r
401     AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
402     TDirectory *savedir = gDirectory;\r
403     fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
404     savedir->cd();\r
405   }\r
406 }\r
407 \r
408 //____________________________________________________________________\r
409 void AliTRDrecoTask::Terminate(Option_t *)\r
410 {\r
411   //\r
412   // Terminate\r
413   //\r
414 \r
415   if(fgDebugStream){ \r
416     delete fgDebugStream;\r
417     fgDebugStream = NULL;\r
418   }\r
419   fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
420   if(fContainer && fRunTerminate){\r
421     PostProcess();\r
422     MakeSummary();\r
423   }\r
424 }\r
425 \r
426 //________________________________________________________\r
427 void AliTRDrecoTask::Adjust(TF1 *f, TH1 * const h)\r
428 {\r
429 // Helper function to avoid duplication of code\r
430 // Make first guesses on the fit parameters\r
431 \r
432   // find the intial parameters of the fit !! (thanks George)\r
433   Int_t nbinsy = Int_t(.5*h->GetNbinsX());\r
434   Double_t sum = 0.;\r
435   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;\r
436   f->SetParLimits(0, 0., 3.*sum);\r
437   f->SetParameter(0, .9*sum);\r
438 \r
439   f->SetParLimits(1, -.2, .2);\r
440   f->SetParameter(1, -0.1);\r
441 \r
442   f->SetParLimits(2, 0., 4.e-1);\r
443   f->SetParameter(2, 2.e-2);\r
444   if(f->GetNpar() <= 4) return;\r
445 \r
446   f->SetParLimits(3, 0., sum);\r
447   f->SetParameter(3, .1*sum);\r
448 \r
449   f->SetParLimits(4, -.3, .3);\r
450   f->SetParameter(4, 0.);\r
451 \r
452   f->SetParLimits(5, 0., 1.e2);\r
453   f->SetParameter(5, 2.e-1);\r
454 }\r