Bug fixed to avoid an infinite loop when no proper TH2 histogram found.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPreprocessor.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 // PHOS Preprocessor class. It runs by Shuttle at the end of the run,
20 // calculates calibration coefficients and dead/bad channels
21 // to be posted in OCDB
22 //
23 // Author: Boris Polichtchouk, 4 October 2006
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include "AliPHOSPreprocessor.h"
27 #include "AliLog.h"
28 #include "AliCDBMetaData.h"
29 #include "AliPHOSEmcCalibData.h"
30 #include "TFile.h"
31 #include "TH1.h"
32 #include "TF1.h"
33 #include "TH2.h"
34 #include "TMap.h"
35 #include "TRandom.h"
36 #include "TKey.h"
37 #include "TList.h"
38 #include "TObjString.h"
39 #include "TObject.h"
40 #include "TString.h"
41 #include "TMath.h"
42 #include "AliPHOSEmcBadChannelsMap.h"
43
44 ClassImp(AliPHOSPreprocessor)
45
46   Double_t rgaus(Double_t *x, Double_t *par)
47 {
48   Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
49                                        (2*par[2]*par[2]) );
50   return gaus;
51 }
52
53 //_______________________________________________________________________________________
54 AliPHOSPreprocessor::AliPHOSPreprocessor() :
55 AliPreprocessor("PHS",0)
56 {
57   //default constructor
58 }
59
60 //_______________________________________________________________________________________
61 AliPHOSPreprocessor::AliPHOSPreprocessor(AliShuttleInterface* shuttle):
62 AliPreprocessor("PHS",shuttle)
63 {
64   // Constructor
65 }
66
67 //_______________________________________________________________________________________
68 UInt_t AliPHOSPreprocessor::Process(TMap* /*valueSet*/)
69 {
70   // process data retrieved by the Shuttle
71   
72   TString runType = GetRunType();
73   Log(Form("Run type: %s",runType.Data()));
74
75   if(runType=="STANDALONE") {
76     Bool_t ledOK = ProcessLEDRun();
77     if(ledOK) return 0;
78     else
79       return 1;
80   }
81   
82   if(runType=="PHYSICS") {
83     
84     Bool_t badmap_OK = FindBadChannelsEmc();
85     if(!badmap_OK) Log(Form("WARNING!! FindBadChannels() completed with BAD status!"));
86     
87     Bool_t calibEmc_OK = CalibrateEmc();
88
89     if(calibEmc_OK && badmap_OK) return 0;
90     else
91       return 1;
92   }
93
94   Log(Form("Unknown run type %s. Do nothing and return OK.",runType.Data()));
95   return 0;
96
97 }
98
99
100 Bool_t AliPHOSPreprocessor::ProcessLEDRun()
101 {
102   //Process LED run, fill bad channels map.
103
104   AliPHOSEmcBadChannelsMap badMap;
105
106   TList* list = GetFileSources(kDAQ, "LED");
107   if(!list) {
108     Log("Sources list for LED run not found, exit.");
109     return kFALSE;
110   }
111
112   TIter iter(list);
113   TObjString *source;
114   char hnam[80];
115   TH1F* histo=0;
116   
117   while ((source = dynamic_cast<TObjString *> (iter.Next()))) {
118
119     AliInfo(Form("found source %s", source->String().Data()));
120
121     TString fileName = GetFile(kDAQ, "LED", source->GetName());
122     AliInfo(Form("Got filename: %s",fileName.Data()));
123
124     TFile f(fileName);
125
126     if(!f.IsOpen()) {
127       Log(Form("File %s is not opened, something goes wrong!",fileName.Data()));
128       return kFALSE;
129     }
130
131     const Int_t nMod=5; // 1:5 modules
132     const Int_t nCol=56; //1:56 columns in each module
133     const Int_t nRow=64; //1:64 rows in each module
134
135     // Check for dead channels    
136     Log(Form("Begin check for dead channels."));
137
138     for(Int_t mod=0; mod<nMod; mod++) {
139       for(Int_t col=0; col<nCol; col++) {
140         for(Int_t row=0; row<nRow; row++) {
141           sprintf(hnam,"mod%dcol%drow%d",mod,col,row);
142           histo = (TH1F*)f.Get(hnam);
143           if(histo)
144             if (histo->GetMean()<1) {
145               Log(Form("Channel: [%d,%d,%d] seems dead, <E>=%.1f.",mod,col,row,histo->GetMean()));
146               badMap.SetBadChannel(mod,col,row);
147             }
148         }
149       }
150     }
151
152   }
153
154   //Store bad channels map
155   AliCDBMetaData badMapMetaData;
156
157   //Bad channels data valid from current run fRun until updated (validityInfinite=kTRUE)
158   Bool_t result = Store("Calib", "EmcBadChannels", &badMap, &badMapMetaData, fRun, kTRUE);
159   return result;
160
161 }
162
163 Float_t AliPHOSPreprocessor::HG2LG(Int_t mod, Int_t X, Int_t Z, TFile* f)
164 {
165   //Calculates High gain to Low gain ratio 
166   //for crystal at the position (X,Z) in the PHOS module mod.
167   
168   char hname[128];
169   sprintf(hname,"%d_%d_%d",mod,X,Z);
170
171   TH1F* h1 = (TH1F*)f->Get(hname);
172   if(!h1) return 16.;
173
174   if(!h1->GetEntries()) return 16.;
175   if(h1->GetMaximum()<10.) return 16.;
176
177   Double_t max = h1->GetBinCenter(h1->GetMaximumBin()); // peak
178   Double_t xmin = max - (h1->GetRMS()/3);
179   Double_t xmax = max + (h1->GetRMS()/2);
180   //       Double_t xmin = max - (h1->GetRMS());
181   //       Double_t xmax = max + (h1->GetRMS());
182
183   TF1* gaus1 = new TF1("gaus1",rgaus,xmin,xmax,3);
184   gaus1->SetParNames("Constant","Mean","Sigma");
185   gaus1->SetParameter("Constant",h1->GetMaximum());
186   gaus1->SetParameter("Mean",max);
187   gaus1->SetParameter("Sigma",1.);
188   gaus1->SetLineColor(kBlue);
189   h1->Fit(gaus1,"LERQ+");
190
191   Log(Form("\t%s: %d entries, mean=%.3f, peak=%.3f, rms= %.3f. HG/LG = %.3f\n",
192            h1->GetTitle(),h1->GetEntries(),h1->GetMean(),max,h1->GetRMS(),
193            gaus1->GetParameter("Mean"))); 
194
195   return gaus1->GetParameter("Mean");
196
197 }
198
199 Bool_t AliPHOSPreprocessor::FindBadChannelsEmc()
200 {
201   //Creates the bad channels map for PHOS EMC.
202
203   // The file fileName contains histograms which have been produced by DA2 detector algorithm.
204   // It is a responsibility of the SHUTTLE framework to form the fileName.
205
206   AliPHOSEmcBadChannelsMap badMap;
207
208   TList* list = GetFileSources(kDAQ, "BAD_CHANNELS");
209
210   if(!list) {
211     Log("Sources list for BAD_CHANNELS not found, exit.");
212     return kFALSE;
213   }
214
215   if(!list->GetEntries()) {
216     Log(Form("Got empty sources list. It seems DA2 did not produce any files!"));
217     return kFALSE;
218   }
219
220   TIter iter(list);
221   TObjString *source;
222   char hnam[80];
223   TH1F* h1=0;
224
225   const Float_t fQualityCut = 1.;
226   Int_t nGoods[5] = {0,0,0,0,0};
227   
228   while ((source = dynamic_cast<TObjString *> (iter.Next()))) {
229
230     AliInfo(Form("found source %s", source->String().Data()));
231
232     TString fileName = GetFile(kDAQ, "BAD_CHANNELS", source->GetName());
233     AliInfo(Form("Got filename: %s",fileName.Data()));
234
235     TFile f(fileName);
236
237     if(!f.IsOpen()) {
238       Log(Form("File %s is not opened, something goes wrong!",fileName.Data()));
239       return kFALSE;
240     }
241     
242     Log(Form("Begin check for bad channels."));
243     
244     for(Int_t mod=0; mod<5; mod++) {
245       for(Int_t iX=0; iX<64; iX++) {
246         for(Int_t iZ=0; iZ<56; iZ++) {
247
248           sprintf(hnam,"%d_%d_%d_%d",mod,iX,iZ,1); // high gain 
249           h1 = (TH1F*)f.Get(hnam);
250
251           if(h1) {
252             Double_t mean = h1->GetMean();
253             
254             if(mean)
255               Log(Form("iX=%d iZ=%d gain=%d   mean=%.3f\n",iX,iZ,1,mean));
256
257             if( mean>0 && mean<fQualityCut ) { 
258               nGoods[mod]++; 
259             }
260             else
261               badMap.SetBadChannel(mod+1,iZ+1,iX+1); //module, col,row
262           }
263
264         }
265       }
266
267       if(nGoods[mod])
268         Log(Form("Module %d: %d good channels.",mod,nGoods[mod]));
269     }
270
271     
272   } // end of loop over sources
273   
274   // Store the bad channels map.
275   
276   AliCDBMetaData md;
277   md.SetResponsible("Boris Polishchuk");
278
279   // Data valid from current run fRun until being updated (validityInfinite=kTRUE)
280   Bool_t result = Store("Calib", "EmcBadChannels", &badMap, &md, fRun, kTRUE);
281   return result;
282
283 }
284
285 Bool_t AliPHOSPreprocessor::CalibrateEmc()
286 {
287   // I.  Calculates the set of calibration coefficients to equalyze the mean energies deposited at high gain.
288   // II. Extracts High_Gain/Low_Gain ratio for each channel.
289
290   // The file fileName contains histograms which have been produced by DA1 detector algorithm.
291   // It is a responsibility of the SHUTTLE framework to form the fileName.
292
293   gRandom->SetSeed(0); //the seed is set to the current  machine clock!
294   AliPHOSEmcCalibData calibData;
295   
296   TList* list = GetFileSources(kDAQ, "AMPLITUDES");
297   
298   if(!list) {
299     Log("Sources list not found, exit.");
300     return 1;
301   }
302
303   if(!list->GetEntries()) {
304     Log(Form("Got empty sources list. It seems DA1 did not produce any files!"));
305     return kFALSE;
306   }
307
308   TIter iter(list);
309   TObjString *source;
310   
311   while ((source = dynamic_cast<TObjString *> (iter.Next()))) {
312     AliInfo(Form("found source %s", source->String().Data()));
313
314     TString fileName = GetFile(kDAQ, "AMPLITUDES", source->GetName());
315     AliInfo(Form("Got filename: %s",fileName.Data()));
316
317     TFile f(fileName);
318
319     if(!f.IsOpen()) {
320       Log(Form("File %s is not opened, something goes wrong!",fileName.Data()));
321       return 1;
322     }
323     
324     const Int_t nMod=5; // 1:5 modules
325     const Int_t nCol=56; //1:56 columns in each module
326     const Int_t nRow=64; //1:64 rows in each module
327
328     Double_t coeff;
329     char hnam[80];
330     TH2F* h2=0;
331     TH1D* h1=0;
332     
333     //Get the reference histogram
334     //(author: Gustavo Conesa Balbastre)
335
336     TList * keylist = f.GetListOfKeys();
337     Int_t nkeys   = f.GetNkeys();
338     Bool_t ok = kFALSE;
339     TKey  *key;
340     Int_t ikey = 0;
341     Int_t counter = 0;
342     TH1D* hRef = 0;
343
344     //Check if the file contains any histogram
345     
346     if(nkeys< 2){
347       Log(Form("Not enough histograms (%d) for calibration.",nkeys));
348       return 1;
349     }
350     
351     while(!ok){
352       ikey = gRandom->Integer(nkeys);
353       key = (TKey*)keylist->At(ikey);
354       TObject* obj = f.Get(key->GetName());
355       TString cname(obj->ClassName());
356       if(cname == "TH2F") {
357         h2 = (TH2F*)obj;
358         TString htitl = h2->GetTitle();
359         if(htitl.Contains("and gain 1")) {
360           hRef = h2->ProjectionX();
361           hRef->SetBins(1000,0.,1000.); // to cut off saturation peak at 1023
362           // Check if the reference histogram has too little statistics
363           if(hRef->GetEntries()>2) ok=kTRUE;
364         }
365       }
366       
367       counter++;
368       
369       if(!ok && counter > nkeys){
370         Log("No histogram with enough statistics for reference. Exit.");
371         return 1;
372       }
373     }
374     
375     Log(Form("reference histogram %s, %.1f entries, mean=%.3f, rms=%.3f.",
376              hRef->GetName(),hRef->GetEntries(),
377              hRef->GetMean(),hRef->GetRMS()));
378
379     Double_t refMean=hRef->GetMean();
380     
381     // Calculates relative calibration coefficients for all non-zero channels
382     
383     for(Int_t mod=0; mod<nMod; mod++) {
384       for(Int_t col=0; col<nCol; col++) {
385         for(Int_t row=0; row<nRow; row++) {
386
387           //High Gain to Low Gain ratio
388           Float_t ratio = HG2LG(mod,row,col,&f);
389           calibData.SetHighLowRatioEmc(mod+1,col+1,row+1,ratio);
390           
391           sprintf(hnam,"%d_%d_%d_1",mod,row,col); // high gain!
392           h2 = (TH2F*)f.Get(hnam);
393
394           //TODO: dead channels exclusion!
395           if(h2) {
396             h1 = h2->ProjectionX();
397             h1->SetBins(1000,0.,1000.); // to cut off saturation peak at 1023
398             coeff = h1->GetMean()/refMean;
399             if(coeff>0)
400               calibData.SetADCchannelEmc(mod+1,col+1,row+1,1./coeff);
401             else 
402               calibData.SetADCchannelEmc(mod+1,col+1,row+1,1.);
403             AliInfo(Form("mod %d col %d row %d  coeff %f\n",mod,col,row,coeff));
404           }
405           else
406             calibData.SetADCchannelEmc(mod+1,col+1,row+1,1.); 
407         }
408       }
409     }
410     
411     f.Close();
412   }
413   
414   //Store EMC calibration data
415   
416   AliCDBMetaData emcMetaData;
417   Bool_t result = Store("Calib", "EmcGainPedestals", &calibData, &emcMetaData, 0, kTRUE); // valid from 0 to infinity);
418   return result;
419
420 }
421
422      
423