.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / TPC / TPCcalib / AliTPCCalibKr.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
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 "AliTPCCalROC.h"
34 #include "AliTPCCalPad.h"
35 #include "AliTPCROC.h"
36 #include "AliMathBase.h"
37 #include "TTreeStream.h"
38
39 //date
40 #include "event.h"
41
42 //header file
43 #include "AliTPCCalibKr.h"
44
45 //----------------------------------------------------------------------------
46 // The AliTPCCalibKr class description (TPC Kr calibration).
47 //
48 //
49 // The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
50 // its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).   
51 // 
52 // The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
53 // These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration 
54 // parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
55 // In addition the debugCalibKr.root file with debug information is created.
56 //
57
58 /*
59  
60 // Usage example:
61 //
62
63 // 1. Analyse output histograms:
64 TFile f("outHistFile.root");
65 AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr");
66 obj->SetRadius(0,0);
67 obj->SetStep(1,1);
68 obj->Analyse();
69
70 // 2. See calibration parameters e.g.:
71 TFile f("calibKr.root");
72 spectrMean->GetCalROC(70)->GetValue(40,40);
73 fitMean->GetCalROC(70)->GetValue(40,40);
74
75 // 3. See debug information e.g.:
76 TFile f("debugCalibKr.root");
77 .ls;
78
79 // -- Print calibKr TTree content 
80 calibKr->Print();
81
82 // -- Draw calibKr TTree variables
83 calibKr.Draw("fitMean");
84
85 */
86
87
88 //
89 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
90 //-----------------------------------------------------------------------------
91
92 ClassImp(AliTPCCalibKr)
93
94 AliTPCCalibKr::AliTPCCalibKr() : 
95   TObject(),
96   fASide(kTRUE),
97   fCSide(kTRUE),
98   fHistoKrArray(72),
99   fADCOverClustSizeMin(0.0),
100   fADCOverClustSizeMax(1.0e9),
101   fMaxADCOverClustADCMin(0.0),
102   fMaxADCOverClustADCMax(1.0e9),
103   fTimeMin(0.0),
104   fTimeMax(1.0e9),
105   fClustSizeMin(0.0),
106   fClustSizeMax(1.0e9),
107   fTimebinRmsIrocMin(0.0),
108   fPadRmsIrocMin(0.0),
109   fRowRmsIrocMin(0.0),
110   fClusterPadSize1DIrocMax(200),
111   fCurveCoefficientIroc(1.0e9),
112   fTimebinRmsOrocMin(0.0),
113   fPadRmsOrocMin(0.0),
114   fRowRmsOrocMin(0.0),
115   fClusterPadSize1DOrocMax(200),
116   fCurveCoefficientOroc(1.0e9),
117   fIrocHistogramMin(100.),
118   fIrocHistogramMax(6000.),
119   fIrocHistogramNbins(200),
120   fOrocHistogramMin(100.),
121   fOrocHistogramMax(5500.),
122   fOrocHistogramNbins(200),
123   fRowRadius(0),
124   fPadRadius(0),
125   fRowStep(1),
126   fPadStep(1)
127
128 {
129   //
130   // default constructor
131   //
132
133   // TObjArray with histograms
134   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
135   fHistoKrArray.Clear();
136  
137   // init cuts (by Stefan)
138 //  SetADCOverClustSizeRange(7,200);
139 //  SetMaxADCOverClustADCRange(0.01,0.4);
140 //  SetTimeRange(200,600);
141 //  SetClustSizeRange(6,200);
142
143   //init  cuts (by Adam)
144   SetTimebinRmsMin(1.6,0.8);
145   SetPadRmsMin(0.825,0.55);
146   SetRowRmsMin(0.1,0.1);
147   SetClusterPadSize1DMax(15,11);
148   SetCurveCoefficient(1.,2.);
149  
150   //set histograms settings
151   SetIrocHistogram(200,100,6000);
152   SetOrocHistogram(200,100,5500);
153
154   // init histograms
155   //Init();
156 }
157
158 //_____________________________________________________________________
159 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : 
160   TObject(pad),
161   
162   fASide(pad.fASide),
163   fCSide(pad.fCSide),
164   fHistoKrArray(72),
165   fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
166   fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
167   fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
168   fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
169   fTimeMin(pad.fTimeMin),
170   fTimeMax(pad.fTimeMax),
171   fClustSizeMin(pad.fClustSizeMin),
172   fClustSizeMax(pad.fClustSizeMax),
173   fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),
174   fPadRmsIrocMin(pad.fPadRmsIrocMin),
175   fRowRmsIrocMin(pad.fRowRmsIrocMin),
176   fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),
177   fCurveCoefficientIroc(pad.fCurveCoefficientIroc),
178   fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),
179   fPadRmsOrocMin(pad.fPadRmsOrocMin),
180   fRowRmsOrocMin(pad.fRowRmsOrocMin),
181   fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),
182   fCurveCoefficientOroc(pad.fCurveCoefficientOroc),
183   fIrocHistogramMin(pad.fIrocHistogramMin),
184   fIrocHistogramMax(pad.fIrocHistogramMax),
185   fIrocHistogramNbins(pad.fIrocHistogramNbins),
186   fOrocHistogramMin(pad.fOrocHistogramMin),
187   fOrocHistogramMax(pad.fOrocHistogramMax),
188   fOrocHistogramNbins(pad.fOrocHistogramNbins),
189   fRowRadius(pad.fRowRadius),
190   fPadRadius(pad.fPadRadius),
191   fRowStep(pad.fRowStep),
192   fPadStep(pad.fPadStep)
193
194 {
195   // copy constructor
196  
197   // TObjArray with histograms
198   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
199   fHistoKrArray.Clear();
200
201   for (Int_t iSec = 0; iSec < 72; ++iSec) 
202   {
203     TH3F *hOld = pad.GetHistoKr(iSec);
204         if(hOld) {
205       TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); 
206       fHistoKrArray.AddAt(hNew,iSec);
207         }
208   }
209 }
210
211 //_____________________________________________________________________
212 AliTPCCalibKr::~AliTPCCalibKr() 
213 {
214   //
215   // destructor
216   //
217
218   // for (Int_t iSec = 0; iSec < 72; ++iSec) {
219   //     if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
220   // }
221   fHistoKrArray.Delete();
222 }
223
224 //_____________________________________________________________________
225 AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
226 {
227   // assignment operator
228
229   if (&source == this) return *this;
230   new (this) AliTPCCalibKr(source);
231
232   return *this;
233 }
234
235 //_____________________________________________________________________
236 void AliTPCCalibKr::Init()
237 {
238   // 
239   // init output histograms 
240   //
241
242   // add histograms to the TObjArray
243   for(Int_t i=0; i<72; ++i) {
244     
245     // C - side
246     if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
247       TH3F *hist = CreateHisto(i);
248       if(hist) fHistoKrArray.AddAt(hist,i);
249     }
250     
251     // A - side
252     if(IsCSide(i) == kFALSE && fASide == kTRUE) {
253       TH3F *hist = CreateHisto(i);
254       if(hist) fHistoKrArray.AddAt(hist,i);
255     }
256   }
257 }
258  
259 //_____________________________________________________________________
260 Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
261 {
262   //
263   // process events 
264   // call event by event
265   //
266
267   if(cluster) Update(cluster);
268   else return kFALSE;
269  
270  return kTRUE;
271 }
272
273 //_____________________________________________________________________
274 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
275 {
276     //
277     // create new histogram
278         //
279     char name[256];
280         TH3F *h;
281
282         snprintf(name,256,"ADCcluster_ch%d",chamber);
283
284     if( IsIROC(chamber) == kTRUE ) 
285         {
286           h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);
287         } else {
288           h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);
289         }
290     h->SetXTitle("padrow");
291     h->SetYTitle("pad");
292     h->SetZTitle("fADC");
293
294 return h;
295 }
296
297 //_____________________________________________________________________
298 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
299 {
300 // check if IROCs
301 // returns kTRUE if IROCs and kFALSE if OROCs 
302
303   if(chamber>=0 && chamber<36) return kTRUE;
304
305 return kFALSE;
306 }
307
308 //_____________________________________________________________________
309 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
310 {
311 // check if C side
312 // returns kTRUE if C side and kFALSE if A side
313
314   if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
315
316 return kFALSE;
317 }
318  
319 //_____________________________________________________________________
320 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
321 {
322   //
323   // fill existing histograms
324   //
325
326   if (!Accept(cl)) return kFALSE;
327   TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
328   if(!h) return kFALSE;
329   
330   h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
331   
332   return kTRUE;
333 }
334
335 //_____________________________________________________________________
336 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
337   //
338   // cuts
339   //
340   /*
341     TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings - 
342     TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
343     TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
344     TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
345     TCut cutR4("cutR4","fMax.fTime>200");   // noise removal
346     TCut cutR5("cutR5","fMax.fTime<600");   // noise removal
347     TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
348     TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
349   */
350   /*
351   //R0
352   if ((float)cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
353   // R1
354   if ((float)cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
355   //R2
356   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
357   //R3
358   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
359   //R4
360   if (cl->GetMax().GetTime() < 200) return kFALSE;
361   //R5
362   if (cl->GetMax().GetTime() > 600) return kFALSE;
363   //S1
364   if (cl->GetSize()>200) return kFALSE;
365   if (cl->GetSize()<6)  return kFALSE;
366
367   SetADCOverClustSizeRange(7,200);
368   SetMaxADCOverClustADCRange(0.01,0.4);
369   SetTimeRange(200,600);
370   SetClustSizeRange(6,200);
371   */
372
373   //R0
374   if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax)        return kFALSE;
375   // R1
376   if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin)          return kFALSE;
377   //R2
378   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax)  return kFALSE;
379   //R3
380   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
381   //R4
382   if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
383   //R5
384   if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
385   //S1
386   if (cl->GetSize()>fClustSizeMax) return kFALSE;
387   if (cl->GetSize()<fClustSizeMin)  return kFALSE;
388
389   //
390   // cuts by Adam
391   //
392   /*
393     TCut cutAI0("cutAI0","fTimebinRMS>1.6");
394     TCut cutAI1("cutAI1","fPadRMS>0.825");  
395     TCut cutAI2("cutAI2","fRowRMS>0.1");    
396     TCut cutAI3("cutAI3","fPads1D<15");  
397     TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0"); 
398
399     TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;
400
401     TCut cutAO0("cutAO0","fTimebinRMS>0.8");
402     TCut cutAO1("cutAO1","fPadRMS>0.55");  
403     TCut cutAO2("cutAO2","fRowRMS>0.1");    
404     TCut cutAO3("cutAO3","fPads1D<11");  
405     TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0"); 
406
407     TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;
408     TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");
409
410   */
411   if(cl->GetSec()<36){  //IROCs
412     //AI0
413     if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;
414     //AI1
415     if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;
416     //AI2
417     if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;
418     //AI3
419     if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;
420     //AI4
421     if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
422   }else{  //OROCs
423     //AO0
424     if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;
425     //AO1
426     if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;
427     //AO2
428     if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;
429     //AO3
430     if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;
431     //AO4
432     if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
433   }
434
435   return kTRUE;
436
437 }
438
439 //_____________________________________________________________________
440 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
441 {
442   // get histograms from fHistoKrArray
443   return (TH3F*) fHistoKrArray.At(chamber);
444 }
445  
446 //_____________________________________________________________________
447 void AliTPCCalibKr::Analyse() 
448 {
449   //
450   // analyse the histograms and extract krypton calibration parameters
451   //
452
453   // AliTPCCalPads that will contain the calibration parameters
454   AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
455   AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
456   AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
457   AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
458   AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
459   AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
460
461   // file stream for debugging purposes
462   TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
463
464   // if entries in spectrum less than minEntries, then the fit won't be performed
465   Int_t minEntries = 1; //300;
466
467   Double_t windowFrac = 0.12;
468   // the 3d histogram will be projected on the pads given by the following window size
469   // set the numbers to 0 if you want to do a pad-by-pad calibration
470   UInt_t rowRadius = fRowRadius;//4;
471   UInt_t padRadius = fPadRadius;//4;
472   
473   // the step size by which pad and row are incremented is given by the following numbers
474   // set them to 1 if you want the finest granularity
475   UInt_t rowStep = fRowStep;//1;//2    // formerly: 2*rowRadius
476   UInt_t padStep = fPadStep;//1;//4    // formerly: 2*padRadius
477
478   for (Int_t chamber = 0; chamber < 72; chamber++) {
479     //if (chamber != 71) continue;
480     AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
481     
482     // Usually I would traverse each pad, take the spectrum for its neighbourhood and
483     // obtain the calibration parameters. This takes very long, so instead I assign the same
484     // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
485     UInt_t nRows = roc.GetNrows();
486     for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
487       UInt_t nPads = roc.GetNPads(iRow);
488       //if (iRow >= 10) break;
489       for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
490         //if (iPad >= 20) break;
491         TH3F* h = GetHistoKr(chamber);
492         if (!h) continue;
493         
494         // the 3d histogram will be projected on the pads given by the following bounds
495         // for rows and pads
496         Int_t rowLow = iRow - rowRadius;
497         UInt_t rowUp = iRow + rowRadius + rowStep-1;
498         Int_t padLow = iPad - padRadius;
499         UInt_t padUp = iPad + padRadius + padStep-1;
500         // if window goes out of chamber
501         if (rowLow < 0) rowLow = 0;
502         if (rowUp >= nRows) rowUp = nRows - 1;
503         if (padLow < 0) padLow = 0;
504         if (padUp >= nPads) padUp = nPads - 1;
505
506         // project the histogram
507         //TH1D* projH = h->ProjectionZ("projH", rowLow+1, rowUp+1, padLow+1, padUp+1); // SLOW
508         TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);
509     
510         // get the number of entries in the spectrum
511         Double_t entries = projH->GetEntries();
512         if (entries < minEntries) { delete projH; continue; }
513         
514         // get the two calibration parameters mean of spectrum and RMS of spectrum
515         Double_t histMean = projH->GetMean();
516         Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
517     
518         // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
519         Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
520                 Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
521                 Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
522         Double_t integCharge =  projH->Integral(minBin,maxBin); 
523
524         Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
525
526         if (fitResult != 0) {
527           Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
528           //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f,  %f\n", chamber, iRow, iPad, entries, maxEntries);
529     
530           delete projH;
531           continue;
532         }
533     
534         // get the two calibration parameters mean of gauss fit and sigma of gauss fit
535         TF1* gausFit = projH->GetFunction("gaus");
536         Double_t fitMean = gausFit->GetParameter(1);
537         Double_t fitRMS = gausFit->GetParameter(2);
538         Int_t numberFitPoints = gausFit->GetNumberFitPoints();
539         if (numberFitPoints == 0) continue;
540         Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
541         delete gausFit;
542         delete projH;
543         if (fitMean <= 0) continue;
544         //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
545     
546         // write the calibration parameters for each pad that the 3d histogram was projected onto
547         // (with considering the step size) to the CalPads
548         // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
549         // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
550         for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
551           if (r < 0 || r >= (Int_t)nRows) continue;
552           UInt_t nPadsR = roc.GetNPads(r);
553           for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
554             if (p < 0 || p >= (Int_t)nPadsR) continue;
555             spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
556             spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
557             fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
558             fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
559             fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
560             entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
561
562             (*debugStream) << "calibKr" <<
563               "sector=" << chamber <<          // chamber number
564               "row=" << r <<                   // row number
565               "pad=" << p <<                   // pad number
566               "histMean=" << histMean <<       // mean of the spectrum
567               "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
568               "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
569               "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
570               "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
571               "entries=" << entries <<         // number of entries for the spectrum
572               "\n";
573           }
574         }
575       }
576     }
577   }
578
579   TFile f("calibKr.root", "recreate");
580   spectrMeanCalPad->Write();
581   spectrRMSCalPad->Write();
582   fitMeanCalPad->Write();
583   fitRMSCalPad->Write();
584   fitNormChi2CalPad->Write();
585   entriesCalPad->Write();
586   f.Close();
587   delete spectrMeanCalPad;
588   delete spectrRMSCalPad;
589   delete fitMeanCalPad;
590   delete fitRMSCalPad;
591   delete fitNormChi2CalPad;
592   delete entriesCalPad;
593   delete debugStream;
594 }
595
596 //_____________________________________________________________________
597 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
598 {
599   // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
600   // replaces TH3F::ProjectZ() to gain more speed
601
602   TAxis* xAxis = histo3D->GetXaxis();
603   TAxis* yAxis = histo3D->GetYaxis();
604   TAxis* zAxis = histo3D->GetZaxis();
605   Double_t zMinVal = zAxis->GetXmin();
606   Double_t zMaxVal = zAxis->GetXmax();
607   
608   Int_t nBinsZ = zAxis->GetNbins();
609   TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
610
611   Int_t nx = xAxis->GetNbins()+2;
612   Int_t ny = yAxis->GetNbins()+2;
613   Int_t bin = 0;
614   Double_t entries = 0.;
615   for (Int_t x = xMin; x <= xMax; x++) {
616     for (Int_t y = yMin; y <= yMax; y++) {
617       for (Int_t z = 0; z <= nBinsZ+1; z++) {
618         bin = x + nx * (y + ny * z);
619         Double_t val = histo3D->GetBinContent(bin);
620         projH->Fill(zAxis->GetBinCenter(z), val);
621         entries += val;
622       }
623     }
624   }
625   projH->SetEntries((Long64_t)entries);
626   return projH;
627 }
628
629 //_____________________________________________________________________
630 Long64_t AliTPCCalibKr::Merge(TCollection* list) {
631 // merge function 
632 //
633
634 if (!list)
635 return 0;
636
637 if (list->IsEmpty())
638 return 1;
639
640 TIterator* iter = list->MakeIterator();
641 TObject* obj = 0;
642
643     iter->Reset();
644     Int_t count=0;
645         while((obj = iter->Next()) != 0)
646         {
647           AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
648           if (entry == 0) continue; 
649
650                 for(int i=0; i<72; ++i) { 
651                   if(IsCSide(i) == kTRUE && fCSide == kTRUE) { 
652                   ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
653                   } 
654
655                   if(IsCSide(i) == kFALSE && fASide == kTRUE) { 
656                     ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
657                   } 
658                 } 
659
660           count++;
661         }
662
663 return count;
664 }