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