High_gain/low_gain calculation added; histogram names changed to satisfy AliPHOSDA1...
[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   // The fileName with the histograms which have been produced by
73   // AliPHOSCalibHistoProducer.
74   // It is a responsibility of the SHUTTLE framework to form the fileName
75   
76   TString runType = GetRunType();
77   Log(Form("Run type: %s",runType.Data()));
78
79   if(runType=="LED") {
80     Bool_t ledOK = ProcessLEDRun();
81     if(ledOK) return 0;
82     else
83       return 1;
84   }
85
86   gRandom->SetSeed(0); //the seed is set to the current  machine clock!
87   AliPHOSEmcCalibData calibData;
88
89   TList* list = GetFileSources(kDAQ, "AMPLITUDES");
90   if(!list) {
91     Log("Sources list not found, exit.");
92     return 1;
93   }
94
95   TIter iter(list);
96   TObjString *source;
97   
98   while ((source = dynamic_cast<TObjString *> (iter.Next()))) {
99     AliInfo(Form("found source %s", source->String().Data()));
100
101     TString fileName = GetFile(kDAQ, "AMPLITUDES", source->GetName());
102     AliInfo(Form("Got filename: %s",fileName.Data()));
103
104     TFile f(fileName);
105
106     if(!f.IsOpen()) {
107       Log(Form("File %s is not opened, something goes wrong!",fileName.Data()));
108       return 1;
109     }
110
111     const Int_t nMod=5; // 1:5 modules
112     const Int_t nCol=56; //1:56 columns in each module
113     const Int_t nRow=64; //1:64 rows in each module
114
115     Double_t coeff;
116     char hnam[80];
117     TH2F* h2=0;
118     TH1D* h1=0;
119     
120     //Get the reference histogram
121     //(author: Gustavo Conesa Balbastre)
122
123     TList * keylist = f.GetListOfKeys();
124     Int_t nkeys   = f.GetNkeys();
125     Bool_t ok = kFALSE;
126     TKey  *key;
127     Int_t ikey = 0;
128     Int_t counter = 0;
129     TH1D* hRef = 0;
130
131     //Check if the file contains any histogram
132     
133     if(nkeys< 2){
134       Log(Form("Not enough histograms (%d) for calibration.",nkeys));
135       return 1;
136     }
137     
138     while(!ok){
139       ikey = gRandom->Integer(nkeys);
140       key = (TKey*)keylist->At(ikey);
141       TObject* obj = f.Get(key->GetName());
142       TString cname(obj->ClassName());
143       if(cname == "TH2F") {
144         h2 = (TH2F*)obj;
145         TString htitl = h2->GetTitle();
146         if(htitl.Contains("and gain 1")) {
147           hRef = h2->ProjectionX();
148           hRef->SetBins(1000,0.,1000.); // to cut off saturation peak at 1023
149           // Check if the reference histogram has too little statistics
150           if(hRef->GetEntries()>2) ok=kTRUE;
151           if(!ok && counter > nkeys){
152             Log("No histogram with enough statistics for reference.");
153             return 1;
154           }
155         }
156       }
157       counter++;
158     }
159     
160     Log(Form("reference histogram %s, %.1f entries, mean=%.3f, rms=%.3f.",
161              hRef->GetName(),hRef->GetEntries(),
162              hRef->GetMean(),hRef->GetRMS()));
163
164     Double_t refMean=hRef->GetMean();
165     
166     // Calculates relative calibration coefficients for all non-zero channels
167     
168     for(Int_t mod=0; mod<nMod; mod++) {
169       for(Int_t col=0; col<nCol; col++) {
170         for(Int_t row=0; row<nRow; row++) {
171
172           //High Gain to Low Gain ratio
173           Float_t ratio = HG2LG(mod,row,col,&f);
174           calibData.SetHighLowRatioEmc(mod+1,col+1,row+1,ratio);
175           
176           sprintf(hnam,"%d_%d_%d_1",mod,row,col); // high gain!
177           h2 = (TH2F*)f.Get(hnam);
178
179           //TODO: dead channels exclusion!
180           if(h2) {
181             h1 = h2->ProjectionX();
182             h1->SetBins(1000,0.,1000.); // to cut off saturation peak at 1023
183             coeff = h1->GetMean()/refMean;
184             if(coeff>0)
185               calibData.SetADCchannelEmc(mod+1,col+1,row+1,1./coeff);
186             else 
187               calibData.SetADCchannelEmc(mod+1,col+1,row+1,1.);
188             AliInfo(Form("mod %d col %d row %d  coeff %f\n",mod,col,row,coeff));
189           }
190           else
191             calibData.SetADCchannelEmc(mod+1,col+1,row+1,1.); 
192         }
193       }
194     }
195     
196     f.Close();
197   }
198   
199   //Store EMC calibration data
200   
201   AliCDBMetaData emcMetaData;
202   Bool_t emcOK = Store("Calib", "EmcGainPedestals", &calibData, &emcMetaData);
203
204   if(emcOK) return 0;
205   else
206     return 1;
207
208 }
209
210
211 Bool_t AliPHOSPreprocessor::ProcessLEDRun()
212 {
213   //Process LED run, fill bad channels map.
214
215   AliPHOSEmcBadChannelsMap badMap;
216
217   TList* list = GetFileSources(kDAQ, "LED");
218   if(!list) {
219     Log("Sources list for LED run not found, exit.");
220     return kFALSE;
221   }
222
223   TIter iter(list);
224   TObjString *source;
225   char hnam[80];
226   TH1F* histo=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, "LED", 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     const Int_t nMod=5; // 1:5 modules
243     const Int_t nCol=56; //1:56 columns in each module
244     const Int_t nRow=64; //1:64 rows in each module
245
246     // Check for dead channels    
247     Log(Form("Begin check for dead channels."));
248
249     for(Int_t mod=0; mod<nMod; mod++) {
250       for(Int_t col=0; col<nCol; col++) {
251         for(Int_t row=0; row<nRow; row++) {
252           sprintf(hnam,"mod%dcol%drow%d",mod,col,row);
253           histo = (TH1F*)f.Get(hnam);
254           if(histo)
255             if (histo->GetMean()<1) {
256               Log(Form("Channel: [%d,%d,%d] seems dead, <E>=%.1f.",mod,col,row,histo->GetMean()));
257               badMap.SetBadChannel(mod,col,row);
258             }
259         }
260       }
261     }
262
263   }
264
265   //Store bad channels map
266   AliCDBMetaData badMapMetaData;
267
268   //Bad channels data valid from current run fRun until updated (validityInfinite=kTRUE)
269   Bool_t result = Store("Calib", "EmcBadChannels", &badMap, &badMapMetaData, fRun, kTRUE);
270   return result;
271
272 }
273
274 Float_t AliPHOSPreprocessor::HG2LG(Int_t mod, Int_t X, Int_t Z, TFile* f)
275 {
276   //Calculates High gain to Low gain ratio 
277   //for crystal at the position (X,Z) in the PHOS module mod.
278   
279   char hname[128];
280   sprintf(hname,"%d_%d_%d",mod,X,Z);
281
282   TH1F* h1 = (TH1F*)f->Get(hname);
283   if(!h1) return 16.;
284
285   if(!h1->GetEntries()) return 16.;
286   if(h1->GetMaximum()<10.) return 16.;
287
288   Double_t max = h1->GetBinCenter(h1->GetMaximumBin()); // peak
289   Double_t xmin = max - (h1->GetRMS()/3);
290   Double_t xmax = max + (h1->GetRMS()/2);
291   //       Double_t xmin = max - (h1->GetRMS());
292   //       Double_t xmax = max + (h1->GetRMS());
293
294   TF1* gaus1 = new TF1("gaus1",rgaus,xmin,xmax,3);
295   gaus1->SetParNames("Constant","Mean","Sigma");
296   gaus1->SetParameter("Constant",h1->GetMaximum());
297   gaus1->SetParameter("Mean",max);
298   gaus1->SetParameter("Sigma",1.);
299   gaus1->SetLineColor(kBlue);
300   h1->Fit(gaus1,"LERQ+");
301
302   Log(Form("\t%s: %d entries, mean=%.3f, peak=%.3f, rms= %.3f. HG/LG = %.3f\n",
303            h1->GetTitle(),h1->GetEntries(),h1->GetMean(),max,h1->GetRMS(),
304            gaus1->GetParameter("Mean"))); 
305
306   return gaus1->GetParameter("Mean");
307
308 }