Bug fix - integer rounding (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibKr.cxx
CommitLineData
72e010d3 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
17//Root includes
18#include <TH1F.h>
19#include <TH1D.h>
20#include <TH2F.h>
21#include <TH3F.h>
22#include <TString.h>
23#include <TMath.h>
24#include <TF1.h>
25#include <TRandom.h>
26#include <TDirectory.h>
27#include <TFile.h>
28#include <TAxis.h>
29//AliRoot includes
30#include "AliRawReader.h"
31#include "AliRawReaderRoot.h"
32#include "AliRawReaderDate.h"
33#include "AliTPCRawStream.h"
34#include "AliTPCCalROC.h"
35#include "AliTPCCalPad.h"
36#include "AliTPCROC.h"
37#include "AliMathBase.h"
38#include "TTreeStream.h"
39#include "AliTPCRawStreamFast.h"
40
41//date
42#include "event.h"
43
44//header file
45#include "AliTPCCalibKr.h"
46
47//----------------------------------------------------------------------------
48// The AliTPCCalibKr class description (TPC Kr calibration).
49//
50//
afb50fd6 51// The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
52// its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).
72e010d3 53//
72e010d3 54// The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
55// These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration
56// parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
57// In addition the debugCalibKr.root file with debug information is created.
58//
afb50fd6 59
60/*
61
72e010d3 62// Usage example:
63//
afb50fd6 64
65// 1. Analyse output histograms:
66TFile f("outHistFile.root");
67AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr")
68obj->Analyse();
69
70// 2. See calibration parameters e.g.:
71TFile f("calibKr.root");
72spectrMean->GetCalROC(70)->GetValue(40,40);
73fitMean->GetCalROC(70)->GetValue(40,40);
74
75// 3. See debug information e.g.:
76TFile f("debugCalibKr.root");
77.ls;
78
72e010d3 79// -- Print calibKr TTree content
afb50fd6 80calibKr->Print();
81
72e010d3 82// -- Draw calibKr TTree variables
afb50fd6 83calibKr.Draw("fitMean");
84
85*/
86
87
72e010d3 88//
89// Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
90//-----------------------------------------------------------------------------
91
92ClassImp(AliTPCCalibKr)
93
94AliTPCCalibKr::AliTPCCalibKr() :
95 TObject(),
96
72e010d3 97 bASide(kTRUE),
98 bCSide(kTRUE),
72e010d3 99 fHistoKrArray(72)
100{
101 //
102 // default constructor
103 //
afb50fd6 104
105 // init histograms
106 Init();
72e010d3 107}
108
109//_____________________________________________________________________
110AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
111 TObject(pad),
112
72e010d3 113 bASide(pad.bASide),
114 bCSide(pad.bCSide),
72e010d3 115 fHistoKrArray(72)
116{
117 // copy constructor
118
119 for (Int_t iSec = 0; iSec < 72; ++iSec)
120 {
121 TH3F *hOld = pad.GetHistoKr(iSec);
122 if(hOld) {
123 TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
124 fHistoKrArray.AddAt(hNew,iSec);
125 }
126 }
127}
128
129//_____________________________________________________________________
130AliTPCCalibKr::~AliTPCCalibKr()
131{
132 //
133 // destructor
134 //
afb50fd6 135
136 for (Int_t iSec = 0; iSec < 72; ++iSec)
137 {
138 TH3F *hNew = (TH3F*)fHistoKrArray.At(iSec);
139 if(hNew) delete hNew; hNew = 0;
140 }
72e010d3 141 fHistoKrArray.Delete();
142}
143
144//_____________________________________________________________________
145AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
146{
147 // assignment operator
148
149 if (&source == this) return *this;
150 new (this) AliTPCCalibKr(source);
151
152 return *this;
153}
154
155//_____________________________________________________________________
156void AliTPCCalibKr::Init()
157{
158 //
159 // init input tree and output histograms
160 //
afb50fd6 161
72e010d3 162 // create output TObjArray
163 fHistoKrArray.Clear();
164
165 // add histograms to the TObjArray
166 for(Int_t i=0; i<72; ++i) {
167
168 // C - side
169 if( IsCSide(i) == kTRUE && bCSide == kTRUE) {
170 TH3F *hist = CreateHisto(i);
171 if(hist) fHistoKrArray.AddAt(hist,i);
172 }
173
174 // A - side
175 if(IsCSide(i) == kFALSE && bASide == kTRUE) {
176 TH3F *hist = CreateHisto(i);
177 if(hist) fHistoKrArray.AddAt(hist,i);
178 }
72e010d3 179 }
180}
72e010d3 181
182//_____________________________________________________________________
afb50fd6 183Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
72e010d3 184{
185 //
186 // process events
187 // call event by event
188 //
189
afb50fd6 190 if(cluster) Update(cluster);
191 else return kFALSE;
72e010d3 192}
193
194//_____________________________________________________________________
195TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
196{
197 //
198 // create new histogram
199 //
200 char name[256];
201 TH3F *h;
202
203 sprintf(name,"ADCcluster_ch%d",chamber);
204
205 if( IsIROC(chamber) == kTRUE )
206 {
207 h = new TH3F(name,name,63,0,63,100,0,100,150,100,3000);
208 } else {
209 h = new TH3F(name,name,96,0,96,100,0,100,150,100,3000);
210 }
211 h->SetXTitle("padrow");
212 h->SetYTitle("pad");
213 h->SetZTitle("fADC");
214
215return h;
216}
217
218//_____________________________________________________________________
219Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
220{
221// check if IROCs
222// returns kTRUE if IROCs and kFALSE if OROCs
223
224 if(chamber>=0 && chamber<36) return kTRUE;
225
226return kFALSE;
227}
228
229//_____________________________________________________________________
230Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
231{
232// check if C side
233// returns kTRUE if C side and kFALSE if A side
234
235 if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
236
237return kFALSE;
238}
239
240//_____________________________________________________________________
241Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
242{
243 //
244 // fill existing histograms
245 //
72e010d3 246
045a55c4 247 if (!Accept(cl)) return kFALSE;
248 TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
249 if(!h) return kFALSE;
250
251 h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
252
253 return kTRUE;
254}
255
afb50fd6 256//_____________________________________________________________________
045a55c4 257Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){
258 //
259 // cuts
260 //
261 /*
262 TCut cutR0("cutR0","fADCcluster/fSize<200"); // adjust it according v seetings -
263 TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
264 TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4"); // digital noise removal
265 TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
266 TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
267 TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
268 */
269 //R0
afb50fd6 270 if (cl->GetADCcluster()/ Float_t(cl->GetSize()) >200) return kFALSE;
045a55c4 271 // R1
afb50fd6 272 if (cl->GetADCcluster()/ Float_t(cl->GetSize()) <7) return kFALSE;
045a55c4 273 //R2
afb50fd6 274 if (cl->GetMax().GetAdc()/ Float_t(cl->GetADCcluster()) >0.4) return kFALSE;
045a55c4 275 //R3
afb50fd6 276 if (cl->GetMax().GetAdc()/ Float_t(cl->GetADCcluster()) <0.01) return kFALSE;
045a55c4 277 //S1
278 if (cl->GetSize()>200) return kFALSE;
279 if (cl->GetSize()<6) return kFALSE;
280 return kTRUE;
72e010d3 281
72e010d3 282}
283
284//_____________________________________________________________________
285TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
286{
287 // get histograms from fHistoKrArray
288 return (TH3F*) fHistoKrArray.At(chamber);
289}
72e010d3 290
291//_____________________________________________________________________
292void AliTPCCalibKr::Analyse()
293{
294 //
295 // analyse the histograms and extract krypton calibration parameters
296 //
297
298 // AliTPCCalPads that will contain the calibration parameters
299 AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
300 AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
301 AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
302 AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
303 AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
304 AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
305
306 // file stream for debugging purposes
307 TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
308
309 // if entries in spectrum less than minEntries, then the fit won't be performed
310 Int_t minEntries = 1; //300;
311
312 Double_t windowFrac = 0.12;
313 // the 3d histogram will be projected on the pads given by the following window size
314 // set the numbers to 0 if you want to do a pad-by-pad calibration
315 UInt_t rowRadius = 5;
316 UInt_t padRadius = 10;
317 // the step size by which pad and row are incremented is given by the following numbers
318 // set them to 1 if you want the finest granularity
319 UInt_t rowStep = 1; // formerly: 2*rowRadius
320 UInt_t padStep = 1; // formerly: 2*padRadius
321
322 for (Int_t chamber = 0; chamber < 72; chamber++) {
323 //if (chamber != 71) continue;
324 AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()
325
326 // Usually I would traverse each pad, take the spectrum for its neighbourhood and
327 // obtain the calibration parameters. This takes very long, so instead I assign the same
328 // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
329 UInt_t nRows = roc.GetNrows();
330 for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
331 UInt_t nPads = roc.GetNPads(iRow);
332 //if (iRow >= 10) break;
333 for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
334 //if (iPad >= 20) break;
335 TH3F* h = GetHistoKr(chamber);
336 if (!h) continue;
337
338 // the 3d histogram will be projected on the pads given by the following bounds
339 // for rows and pads
340 Int_t rowLow = iRow - rowRadius;
341 UInt_t rowUp = iRow + rowRadius;
342 Int_t padLow = iPad - padRadius;
343 UInt_t padUp = iPad + padRadius;
344 // if window goes out of chamber
345 if (rowLow < 0) rowLow = 0;
346 if (rowUp >= nRows) rowUp = nRows - 1;
347 if (padLow < 0) padLow = 0;
348 if (padUp >= nPads) padUp = nPads - 1;
349
350 // project the histogram
351 //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
352 TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
353
354 // get the number of entries in the spectrum
355 Double_t entries = projH->GetEntries();
356 if (entries < minEntries) { delete projH; continue; }
357
358 // get the two calibration parameters mean of spectrum and RMS of spectrum
359 Double_t histMean = projH->GetMean();
360 Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
361
362 // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
363 Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
364 Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
365 if (fitResult != 0) {
366 Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i.", chamber, rowLow, rowUp, padLow, padUp);
367 delete projH;
368 continue;
369 }
370
371 // get the two calibration parameters mean of gauss fit and sigma of gauss fit
372 TF1* gausFit = projH->GetFunction("gaus");
373 Double_t fitMean = gausFit->GetParameter(1);
374 Double_t fitRMS = gausFit->GetParameter(2);
375 Int_t numberFitPoints = gausFit->GetNumberFitPoints();
376 if (numberFitPoints == 0) continue;
377 Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
378 delete gausFit;
379 delete projH;
380 if (fitMean <= 0) continue;
afb50fd6 381 //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
72e010d3 382
383 // write the calibration parameters for each pad that the 3d histogram was projected onto
384 // (with considering the step size) to the CalPads
385 // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
386 // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
387 for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
388 if (r < 0 || r >= (Int_t)nRows) continue;
389 UInt_t nPads = roc.GetNPads(r);
390 for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
391 if (p < 0 || p >= (Int_t)nPads) continue;
392 spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
393 spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
394 fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
395 fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
396 fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
397 entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
398
399 (*debugStream) << "calibKr" <<
400 "sector=" << chamber << // chamber number
401 "row=" << r << // row number
402 "pad=" << p << // pad number
403 "histMean=" << histMean << // mean of the spectrum
404 "histRMS=" << histRMS << // RMS of the spectrum divided by the mean
405 "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak
406 "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak
407 "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit
408 "entries=" << entries << // number of entries for the spectrum
409 "\n";
410 }
411 }
412 }
413 }
414 }
415
416 TFile f("calibKr.root", "recreate");
417 spectrMeanCalPad->Write();
418 spectrRMSCalPad->Write();
419 fitMeanCalPad->Write();
420 fitRMSCalPad->Write();
421 fitNormChi2CalPad->Write();
422 entriesCalPad->Write();
423 f.Close();
424 delete spectrMeanCalPad;
425 delete spectrRMSCalPad;
426 delete fitMeanCalPad;
427 delete fitRMSCalPad;
428 delete fitNormChi2CalPad;
429 delete entriesCalPad;
430 delete debugStream;
431}
432
433//_____________________________________________________________________
434TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
435{
436 // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
437 // replaces TH3F::ProjectZ() to gain more speed
438
439 TAxis* xAxis = histo3D->GetXaxis();
440 TAxis* yAxis = histo3D->GetYaxis();
441 TAxis* zAxis = histo3D->GetZaxis();
442 Double_t zMinVal = zAxis->GetXmin();
443 Double_t zMaxVal = zAxis->GetXmax();
444
445 Int_t nBinsZ = zAxis->GetNbins();
446 TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
447
448 Int_t nx = xAxis->GetNbins()+2;
449 Int_t ny = yAxis->GetNbins()+2;
450 Int_t bin = 0;
451 Double_t entries = 0.;
452 for (Int_t x = xMin; x <= xMax; x++) {
453 for (Int_t y = yMin; y <= yMax; y++) {
454 for (Int_t z = 0; z <= nBinsZ+1; z++) {
455 bin = x + nx * (y + ny * z);
456 Double_t val = histo3D->GetBinContent(bin);
457 projH->Fill(zAxis->GetBinCenter(z), val);
458 entries += val;
459 }
460 }
461 }
462 projH->SetEntries((Long64_t)entries);
463 return projH;
464}
afb50fd6 465
466//_____________________________________________________________________
467Long64_t AliTPCCalibKr::Merge(TCollection* list) {
468// merge function
469//
470cout << "Merge " << endl;
471
472if (!list)
473return 0;
474
475if (list->IsEmpty())
476return 1;
477
478TIterator* iter = list->MakeIterator();
479TObject* obj = 0;
480
481 iter->Reset();
482 Int_t count=0;
483 while((obj = iter->Next()) != 0)
484 {
485 AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
486 if (entry == 0) continue;
487
488 for(int i=0; i<72; ++i) {
489 if(IsCSide(i) == kTRUE && bCSide == kTRUE) {
490 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
491 }
492
493 if(IsCSide(i) == kFALSE && bASide == kTRUE) {
494 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
495 }
496 }
497
498 count++;
499 }
500
501return count;
502}