Bug fix. Removed delete statement
[u/mrichter/AliRoot.git] / ESDCheck / AliFMDQATask.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //_________________________________________________________________________
19 // An analysis task to check the FMD data in simulated data
20 //
21 //*-- Hans Hjersing Dalsgaard 
22 //////////////////////////////////////////////////////////////////////////////
23
24 #include <TCanvas.h> 
25 #include <TChain.h>
26 #include <TF1.h> 
27 #include <TFile.h> 
28 #include <TH1D.h> 
29 #include <TROOT.h>
30 #include <TString.h> 
31
32 #include "AliFMDQATask.h" 
33 #include "AliESD.h" 
34 #include "AliLog.h"
35
36 //______________________________________________________________________________
37 AliFMDQATask::AliFMDQATask(const char *name) : 
38   AliAnalysisTask(name,""),  
39   fChain(0),
40   fESD(0), 
41   fOutputContainer(0),
42   fhFMD1i(0),
43   fhFMD2i(0), 
44   fhFMD2o(0), 
45   fhFMD3i(0), 
46   fhFMD3o(0) 
47 {
48   // Constructor.
49   // Input slot #0 works with an Ntuple
50   DefineInput(0, TChain::Class());
51   // Output slot #0 writes into a TH1 container
52   DefineOutput(0,  TObjArray::Class()) ; 
53 }
54
55 //______________________________________________________________________________
56 AliFMDQATask::~AliFMDQATask()
57 {
58   // dtor
59
60   fOutputContainer->Clear() ; 
61   delete fOutputContainer ; 
62
63   delete fhFMD1i ;
64   delete fhFMD2i ; 
65   delete fhFMD2o ; 
66   delete fhFMD3i ; 
67   delete fhFMD3o ;
68
69 }
70
71 //______________________________________________________________________________
72 void AliFMDQATask::ConnectInputData(const Option_t*)
73 {
74   // Initialisation of branch container and histograms 
75     
76   AliInfo(Form("*** Initialization of %s", GetName())) ; 
77   
78   // Get input data
79   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
80   if (!fChain) {
81     AliError(Form("Input 0 for %s not found\n", GetName()));
82     return ;
83   }
84   
85   // One should first check if the branch address was taken by some other task
86   char ** address = (char **)GetBranchAddress(0, "ESD");
87   if (address) {
88     fESD = (AliESD*)(*address);
89   } else {
90     fESD = new AliESD();
91     SetBranchAddress(0, "ESD", &fESD);
92   }
93 }
94
95 //________________________________________________________________________
96 void AliFMDQATask::CreateOutputObjects()
97 {  
98   // create histograms 
99
100   OpenFile(0) ; 
101
102   fhFMD1i = new TH1D("FMD1i", "FMD1i", 100, -0.5, 3);
103   fhFMD2i = new TH1D("FMD2i", "FMD2i", 100, -0.5, 3);
104   fhFMD2o = new TH1D("FMD2o", "FMD2o", 100, -0.5, 3);
105   fhFMD3i = new TH1D("FMD3i", "FMD3i", 100, -0.5, 3);
106   fhFMD3o = new TH1D("FMD3o", "FMD3o", 100, -0.5, 3);
107
108   // create output container
109   
110   fOutputContainer = new TObjArray(5) ; 
111   fOutputContainer->SetName(GetName()) ; 
112
113   fOutputContainer->AddAt(fhFMD1i,             0) ; 
114   fOutputContainer->AddAt(fhFMD2i,             1) ; 
115   fOutputContainer->AddAt(fhFMD2o,             2) ; 
116   fOutputContainer->AddAt(fhFMD3i,             3) ; 
117   fOutputContainer->AddAt(fhFMD3o,             4) ; 
118 }
119
120 //______________________________________________________________________________
121 void AliFMDQATask::Exec(Option_t *) 
122 {
123   // Processing of one event
124     
125   Long64_t entry = fChain->GetReadEntry() ;
126   
127   if (!fESD) {
128     AliError("fESD is not connected to the input!") ; 
129     return ; 
130   }
131   
132   TFile * currentFile = (dynamic_cast<TChain *>(fChain))->GetFile() ; 
133   if ( !((entry-1)%100) ) 
134     AliInfo(Form("%s ----> Processing event # %lld", currentFile->GetName(), entry)) ; 
135             
136   // ************************  FMD *************************************
137  AliESDFMD * fmd = fESD->GetFMDData() ;
138  if ( ! fmd ) {
139    AliError("No FMD found in ESD") ; 
140    return ; 
141  }
142  fmd->CheckNeedUShort(currentFile);
143   Int_t nFMD1 = 0, nFMD2i = 0, nFMD2o = 0, nFMD3i = 0, nFMD3o = 0 ;
144   
145   UShort_t detector = 1 ;
146   for(detector = 1 ; detector <= fmd->MaxDetectors() ; detector++) {
147     Char_t ring = 'O' ;
148     UShort_t sector ; 
149     for(sector = 0 ;sector < fmd->MaxSectors() ; sector++) {
150       UShort_t strip ;
151       for(strip = 0 ; strip < fmd->MaxStrips(); strip++) {
152         if(fmd->Multiplicity(detector, ring, sector, strip) != AliESDFMD::kInvalidMult)
153           RingSelector(detector, ring, fmd->Multiplicity(detector, ring, sector, strip)) ;
154         if( (fmd->Multiplicity(detector, ring, sector, strip) == AliESDFMD::kInvalidMult) && detector == 2 )
155           nFMD2o++ ;
156         if( (fmd->Multiplicity(detector, ring, sector, strip)==AliESDFMD::kInvalidMult) && detector == 3 )
157           nFMD3o++ ;
158       }
159     }
160     ring='I';
161     for(sector = 0; sector < fmd->MaxSectors(); sector++) {
162       UShort_t strip ;
163       for(strip = 0 ; strip < fmd->MaxStrips() ; strip++) {
164         if(fmd->Multiplicity(detector, ring, sector, strip) != AliESDFMD::kInvalidMult)
165           RingSelector(detector, ring, fmd->Multiplicity(detector, ring, sector, strip));
166         if( (fmd->Multiplicity(detector, ring, sector, strip) == AliESDFMD::kInvalidMult) && detector == 1 )
167           nFMD1++;
168         if( (fmd->Multiplicity(detector, ring, sector, strip) == AliESDFMD::kInvalidMult) && detector == 2 )
169           nFMD2i++;
170         if( (fmd->Multiplicity(detector, ring, sector, strip) == AliESDFMD::kInvalidMult) && detector == 3 )
171           nFMD3i++;
172       }     
173     }
174   }
175
176   if(nFMD1>100+10240)
177     AliWarning(Form("number of missing strips in FMD1i too high in event number %lld in file", entry, fChain->GetCurrentFile()->GetName())) ;
178   if(nFMD2i>100+10240)
179     AliWarning(Form("number of missing strips in FMD2i too high in event number %lld in file", entry, fChain->GetCurrentFile()->GetName())) ;
180   if(nFMD2o>100+10240)
181     AliWarning(Form("number of missing strips in FMD2o too high in event number %lld in file", entry, fChain->GetCurrentFile()->GetName())) ;
182   if(nFMD3i>100+10240)
183     AliWarning(Form("number of missing strips in FMD3i too high in event number %lld in file", entry, fChain->GetCurrentFile()->GetName())) ;
184   if(nFMD3o>100+10240)
185     AliWarning(Form("number of missing strips in FMD3o too high in event number %lld in file", entry, fChain->GetCurrentFile()->GetName())) ;
186   
187   PostData(0, fOutputContainer);
188 }
189
190 //______________________________________________________________________________
191 void AliFMDQATask::Terminate(Option_t *)
192 {
193   // Processing when the event loop is ended
194  
195   fOutputContainer = (TObjArray*)GetOutputData(0);
196   if ( ! fOutputContainer ) 
197     return ; 
198
199   fhFMD1i = (TH1D*)fOutputContainer->At(0);
200   fhFMD2i = (TH1D*)fOutputContainer->At(1);
201   fhFMD2o = (TH1D*)fOutputContainer->At(2);
202   fhFMD3i = (TH1D*)fOutputContainer->At(3);
203   fhFMD3o = (TH1D*)fOutputContainer->At(4);
204
205   TCanvas * cFMD1 = new TCanvas("cFMD1", "FMD ESD Test", 400, 10, 600, 700);
206   cFMD1->Divide(3, 2) ; 
207
208   cFMD1->cd(1) ;; 
209   fhFMD1i->Draw() ; 
210   
211   cFMD1->cd(2) ;; 
212   fhFMD2i->Draw() ; 
213
214   cFMD1->cd(3) ;; 
215   fhFMD2o->Draw() ; 
216
217   cFMD1->cd(4) ;; 
218   fhFMD3i->Draw() ; 
219
220   cFMD1->cd(5) ;; 
221   fhFMD3o->Draw() ; 
222   
223   cFMD1->Print("FMD.eps") ;
224
225   Bool_t rv1i = TestHisto(fhFMD1i) ;
226   Bool_t rv2i = TestHisto(fhFMD2i) ;
227   Bool_t rv2o = TestHisto(fhFMD2o) ;
228   Bool_t rv3i = TestHisto(fhFMD2i) ;
229   Bool_t rv3o = TestHisto(fhFMD2o) ;
230   
231   if (  !(rv1i * rv2i * rv2o * rv3i * rv3o) )
232     AliWarning("Possible problem in file !!! Check output!") ;
233
234
235
236   char line[1024] ; 
237   sprintf(line, ".!tar -zcvf %s.tar.gz *.eps", GetName()) ; 
238   gROOT->ProcessLine(line);
239   sprintf(line, ".!rm -fR *.eps"); 
240   gROOT->ProcessLine(line);
241  
242   AliInfo(Form("!!! All the eps files are in %s.tar.gz !!! \n", GetName())) ;
243 }
244
245 //______________________________________________________________________________
246 void AliFMDQATask::RingSelector(const UShort_t detector, const Char_t ring, const Float_t mult) const 
247
248   // fill the histograms for each ring in each detector layer
249
250   if(ring == 'I' && detector == 1)
251     fhFMD1i->Fill(mult) ;
252   if(ring == 'I' && detector == 2)
253     fhFMD2i->Fill(mult) ;
254   if(ring == 'O' && detector == 2)
255     fhFMD2o ->Fill(mult) ;
256   if(ring == 'I' && detector == 3)
257     fhFMD3i ->Fill(mult) ;
258   if(ring == 'O' && detector == 3)
259     fhFMD3o ->Fill(mult) ;
260 }
261 //______________________________________________________________________________
262 Bool_t AliFMDQATask::TestHisto(TH1D * hTest) const 
263 {  
264   // analyses the histogram with a Landau function
265
266   Float_t chiMax = 3, chiLow=0.5 ;
267   Float_t chiSq ;
268   Float_t mpv ;
269   Int_t   ndf ;
270   
271   FitAll(hTest, chiSq, ndf, mpv, chiMax, chiLow);
272
273   if( (chiSq > chiMax || chiSq < chiLow) || mpv < 0.6 || mpv > 1 ) {
274     hTest->Rebin(2) ;
275     FitAll(hTest, chiSq, ndf, mpv, chiMax, chiLow) ;
276   }
277
278   Bool_t   ret   = kFALSE ;
279   const Char_t * test  = "not OK";
280   const Char_t * test2 = "not OK";
281   
282   if(chiSq < chiMax && chiSq > chiLow)
283     test = "OK" ;
284   if(mpv > 0.6 && mpv < 1)
285     test2 = "OK" ;
286  
287   if(!strcmp(test,"OK") && !strcmp(test2,"OK"))
288     ret = kTRUE;
289   
290   if(!strcmp(test,"not OK") || !strcmp(test2,"not OK")) {
291     AliWarning("Bad fit results") ; 
292     printf("Detector : %s\n", hTest->GetName()) ;
293     printf("Landau fit Chi Square / NDF = %f / %d which is %s\n", chiSq*ndf, ndf, test) ; 
294     printf("Landau fit MPV is: %f which is %s\n", mpv, test2) ;
295   }
296   return ret;  
297 }
298 //______________________________________________________________________________
299 void AliFMDQATask::FitAll(TH1D* hTest, Float_t &chiSq, Int_t &ndf, Float_t &mpv, Float_t chiMax, Float_t chiLow ) const 
300 {
301   // fit a histogram with a Landau distribution and returns chi square
302
303   Float_t fitMax = hTest->GetXaxis()->GetXmax() ;
304
305   TH1D hTmp = *hTest ;
306   hTmp.SetAxisRange(0.4,fitMax) ;
307   Int_t   maxBin = hTmp.GetMaximumBin();
308   Float_t max    = hTmp.GetBinCenter(maxBin);
309  
310   hTest->Fit("landau", "QOI", "", max-0.3, fitMax) ;
311   TF1* fitfunc = hTest->GetFunction("landau") ;
312   chiSq = fitfunc->GetChisquare() / fitfunc->GetNDF() ;
313   mpv   = fitfunc->GetParameter(1) ;
314   ndf   = fitfunc->GetNDF() ;
315
316   if( ( chiSq > chiMax || chiSq < chiLow ) || ( mpv < 0.6 || mpv > 1 ) ) {
317     hTest->Fit("landau", "QOI", "", max-0.2, fitMax) ;
318     fitfunc = hTest->GetFunction("landau") ;
319     chiSq   = fitfunc->GetChisquare() / fitfunc->GetNDF() ;
320     mpv     = fitfunc->GetParameter(1) ;
321     ndf     = fitfunc->GetNDF() ;
322   }
323   if( ( chiSq >chiMax || chiSq < chiLow ) || ( mpv < 0.6 || mpv > 1 ) ) {
324     hTest->Fit("landau", "QOI", "", max-0.1, fitMax) ;
325     fitfunc = hTest->GetFunction("landau") ;
326     chiSq   = fitfunc->GetChisquare() / fitfunc->GetNDF() ;
327     mpv     = fitfunc->GetParameter(1) ;
328     ndf     = fitfunc->GetNDF() ;
329   }
330   if( ( chiSq > chiMax || chiSq <chiLow ) || ( mpv < 0.6 || mpv > 1 ) ) {
331     hTest->Fit("landau", "QOI", "", max, fitMax) ;
332     fitfunc = hTest->GetFunction("landau") ;
333     chiSq   = fitfunc->GetChisquare() / fitfunc->GetNDF() ;
334     mpv     = fitfunc->GetParameter(1) ;
335     ndf     = fitfunc->GetNDF(); 
336   }
337   if( ( chiSq > chiMax || chiSq < chiLow ) || ( mpv < 0.6 ||mpv > 1 ) ) {
338     hTest->Fit("landau", "QOI", "", max+0.1, fitMax) ;
339     fitfunc = hTest->GetFunction("landau") ;
340     chiSq   = fitfunc->GetChisquare() / fitfunc->GetNDF() ;
341     mpv     = fitfunc->GetParameter(1) ;
342     ndf     = fitfunc->GetNDF() ;
343   } 
344 }