Just something left from v2...now clean
[u/mrichter/AliRoot.git] / TPC / 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 "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 //
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()).   
53 // 
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 //
59
60 /*
61  
62 // Usage example:
63 //
64
65 // 1. Analyse output histograms:
66 TFile f("outHistFile.root");
67 AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr")
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 {
108   //
109   // default constructor
110   //
111
112   // TObjArray with histograms
113   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
114   fHistoKrArray.Clear();
115  
116   // init cuts
117   SetADCOverClustSizeRange(7,200);
118   SetMaxADCOverClustADCRange(0.01,0.4);
119   SetTimeRange(200,600);
120   SetClustSizeRange(6,200);
121  
122   // init histograms
123   Init();
124 }
125
126 //_____________________________________________________________________
127 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : 
128   TObject(pad),
129   
130   fASide(pad.fASide),
131   fCSide(pad.fCSide),
132   fHistoKrArray(72),
133   fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
134   fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
135   fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
136   fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
137   fTimeMin(pad.fTimeMin),
138   fTimeMax(pad.fTimeMax),
139   fClustSizeMin(pad.fClustSizeMin),
140   fClustSizeMax(pad.fClustSizeMax)
141 {
142   // copy constructor
143  
144   // TObjArray with histograms
145   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
146   fHistoKrArray.Clear();
147
148   for (Int_t iSec = 0; iSec < 72; ++iSec) 
149   {
150     TH3F *hOld = pad.GetHistoKr(iSec);
151         if(hOld) {
152       TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); 
153       fHistoKrArray.AddAt(hNew,iSec);
154         }
155   }
156 }
157
158 //_____________________________________________________________________
159 AliTPCCalibKr::~AliTPCCalibKr() 
160 {
161   //
162   // destructor
163   //
164
165   // for (Int_t iSec = 0; iSec < 72; ++iSec) {
166   //     if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
167   // }
168   fHistoKrArray.Delete();
169 }
170
171 //_____________________________________________________________________
172 AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
173 {
174   // assignment operator
175
176   if (&source == this) return *this;
177   new (this) AliTPCCalibKr(source);
178
179   return *this;
180 }
181
182 //_____________________________________________________________________
183 void AliTPCCalibKr::Init()
184 {
185   // 
186   // init output histograms 
187   //
188
189   // add histograms to the TObjArray
190   for(Int_t i=0; i<72; ++i) {
191     
192         // C - side
193         if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
194       TH3F *hist = CreateHisto(i);
195       if(hist) fHistoKrArray.AddAt(hist,i);
196         }
197     
198         // A - side
199         if(IsCSide(i) == kFALSE && fASide == kTRUE) {
200       TH3F *hist = CreateHisto(i);
201       if(hist) fHistoKrArray.AddAt(hist,i);
202         }
203   }
204 }
205  
206 //_____________________________________________________________________
207 Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
208 {
209   //
210   // process events 
211   // call event by event
212   //
213
214   if(cluster) Update(cluster);
215   else return kFALSE;
216  
217  return kTRUE;
218 }
219
220 //_____________________________________________________________________
221 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
222 {
223     //
224     // create new histogram
225         //
226     char name[256];
227         TH3F *h;
228
229     sprintf(name,"ADCcluster_ch%d",chamber);
230
231     if( IsIROC(chamber) == kTRUE ) 
232         {
233            h = new TH3F(name,name,63,0,63,108,0,108,200,100,2500);
234         } else {
235            h = new TH3F(name,name,96,0,96,140,0,140,200,100,1700);
236         }
237     h->SetXTitle("padrow");
238     h->SetYTitle("pad");
239     h->SetZTitle("fADC");
240
241 return h;
242 }
243
244 //_____________________________________________________________________
245 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
246 {
247 // check if IROCs
248 // returns kTRUE if IROCs and kFALSE if OROCs 
249
250   if(chamber>=0 && chamber<36) return kTRUE;
251
252 return kFALSE;
253 }
254
255 //_____________________________________________________________________
256 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
257 {
258 // check if C side
259 // returns kTRUE if C side and kFALSE if A side
260
261   if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
262
263 return kFALSE;
264 }
265  
266 //_____________________________________________________________________
267 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
268 {
269   //
270   // fill existing histograms
271   //
272
273   if (!Accept(cl)) return kFALSE;
274   TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
275   if(!h) return kFALSE;
276   
277   h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
278   
279   return kTRUE;
280 }
281
282 //_____________________________________________________________________
283 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
284   //
285   // cuts
286   //
287   /*
288     TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings - 
289     TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
290     TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
291     TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
292     TCut cutR4("cutR4","fMax.fTime>200");   // noise removal
293     TCut cutR5("cutR5","fMax.fTime<600");   // noise removal
294     TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
295     TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
296   */
297   /*
298   //R0
299   if ((float)cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
300   // R1
301   if ((float)cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
302   //R2
303   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
304   //R3
305   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
306   //R4
307   if (cl->GetMax().GetTime() < 200) return kFALSE;
308   //R5
309   if (cl->GetMax().GetTime() > 600) return kFALSE;
310   //S1
311   if (cl->GetSize()>200) return kFALSE;
312   if (cl->GetSize()<6)  return kFALSE;
313
314   SetADCOverClustSizeRange(7,200);
315   SetMaxADCOverClustADCRange(0.01,0.4);
316   SetTimeRange(200,600);
317   SetClustSizeRange(6,200);
318   */
319
320   //R0
321   if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax)        return kFALSE;
322   // R1
323   if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin)          return kFALSE;
324   //R2
325   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax)  return kFALSE;
326   //R3
327   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
328   //R4
329   if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
330   //R5
331   if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
332   //S1
333   if (cl->GetSize()>fClustSizeMax) return kFALSE;
334   if (cl->GetSize()<fClustSizeMin)  return kFALSE;
335
336   return kTRUE;
337
338 }
339
340 //_____________________________________________________________________
341 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
342 {
343   // get histograms from fHistoKrArray
344   return (TH3F*) fHistoKrArray.At(chamber);
345 }
346  
347 //_____________________________________________________________________
348 void AliTPCCalibKr::Analyse() 
349 {
350   //
351   // analyse the histograms and extract krypton calibration parameters
352   //
353
354   // AliTPCCalPads that will contain the calibration parameters
355   AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
356   AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
357   AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
358   AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
359   AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
360   AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
361
362   // file stream for debugging purposes
363   TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
364
365   // if entries in spectrum less than minEntries, then the fit won't be performed
366   Int_t minEntries = 1; //300;
367
368   Double_t windowFrac = 0.12;
369   // the 3d histogram will be projected on the pads given by the following window size
370   // set the numbers to 0 if you want to do a pad-by-pad calibration
371   UInt_t rowRadius = 2;
372   UInt_t padRadius = 4;
373   
374   // the step size by which pad and row are incremented is given by the following numbers
375   // set them to 1 if you want the finest granularity
376   UInt_t rowStep = 1;    // formerly: 2*rowRadius
377   UInt_t padStep = 1;    // formerly: 2*padRadius
378
379   for (Int_t chamber = 0; chamber < 72; chamber++) {
380     //if (chamber != 71) continue;
381     AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
382     
383     // Usually I would traverse each pad, take the spectrum for its neighbourhood and
384     // obtain the calibration parameters. This takes very long, so instead I assign the same
385     // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
386     UInt_t nRows = roc.GetNrows();
387     for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
388       UInt_t nPads = roc.GetNPads(iRow);
389       //if (iRow >= 10) break;
390       for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
391         //if (iPad >= 20) break;
392         TH3F* h = GetHistoKr(chamber);
393         if (!h) continue;
394         
395         // the 3d histogram will be projected on the pads given by the following bounds
396         // for rows and pads
397         Int_t rowLow = iRow - rowRadius;
398         UInt_t rowUp = iRow + rowRadius;
399         Int_t padLow = iPad - padRadius;
400         UInt_t padUp = iPad + padRadius;
401         // if window goes out of chamber
402         if (rowLow < 0) rowLow = 0;
403         if (rowUp >= nRows) rowUp = nRows - 1;
404         if (padLow < 0) padLow = 0;
405         if (padUp >= nPads) padUp = nPads - 1;
406
407         // project the histogram
408         //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
409         TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
410     
411         // get the number of entries in the spectrum
412         Double_t entries = projH->GetEntries();
413         if (entries < minEntries) { delete projH; continue; }
414         
415         // get the two calibration parameters mean of spectrum and RMS of spectrum
416         Double_t histMean = projH->GetMean();
417         Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
418     
419         // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
420         Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
421                 Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
422                 Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
423         Double_t integCharge =  projH->Integral(minBin,maxBin); 
424
425         Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
426
427         if (fitResult != 0) {
428           Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
429           //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f,  %f\n", chamber, iRow, iPad, entries, maxEntries);
430     
431           delete projH;
432           continue;
433         }
434     
435         // get the two calibration parameters mean of gauss fit and sigma of gauss fit
436         TF1* gausFit = projH->GetFunction("gaus");
437         Double_t fitMean = gausFit->GetParameter(1);
438         Double_t fitRMS = gausFit->GetParameter(2);
439         Int_t numberFitPoints = gausFit->GetNumberFitPoints();
440         if (numberFitPoints == 0) continue;
441         Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
442         delete gausFit;
443         delete projH;
444         if (fitMean <= 0) continue;
445         //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
446     
447         // write the calibration parameters for each pad that the 3d histogram was projected onto
448         // (with considering the step size) to the CalPads
449         // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
450         // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
451         for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
452           if (r < 0 || r >= (Int_t)nRows) continue;
453           UInt_t nPads = roc.GetNPads(r);
454           for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
455             if (p < 0 || p >= (Int_t)nPads) continue;
456             spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
457             spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
458             fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
459             fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
460             fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
461             entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
462
463             (*debugStream) << "calibKr" <<
464               "sector=" << chamber <<          // chamber number
465               "row=" << r <<                   // row number
466               "pad=" << p <<                   // pad number
467               "histMean=" << histMean <<       // mean of the spectrum
468               "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
469               "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
470               "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
471               "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
472               "entries=" << entries <<         // number of entries for the spectrum
473               "\n";
474           }
475         }
476       }
477     }
478   }
479
480   TFile f("calibKr.root", "recreate");
481   spectrMeanCalPad->Write();
482   spectrRMSCalPad->Write();
483   fitMeanCalPad->Write();
484   fitRMSCalPad->Write();
485   fitNormChi2CalPad->Write();
486   entriesCalPad->Write();
487   f.Close();
488   delete spectrMeanCalPad;
489   delete spectrRMSCalPad;
490   delete fitMeanCalPad;
491   delete fitRMSCalPad;
492   delete fitNormChi2CalPad;
493   delete entriesCalPad;
494   delete debugStream;
495 }
496
497 //_____________________________________________________________________
498 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
499 {
500   // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
501   // replaces TH3F::ProjectZ() to gain more speed
502
503   TAxis* xAxis = histo3D->GetXaxis();
504   TAxis* yAxis = histo3D->GetYaxis();
505   TAxis* zAxis = histo3D->GetZaxis();
506   Double_t zMinVal = zAxis->GetXmin();
507   Double_t zMaxVal = zAxis->GetXmax();
508   
509   Int_t nBinsZ = zAxis->GetNbins();
510   TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
511
512   Int_t nx = xAxis->GetNbins()+2;
513   Int_t ny = yAxis->GetNbins()+2;
514   Int_t bin = 0;
515   Double_t entries = 0.;
516   for (Int_t x = xMin; x <= xMax; x++) {
517     for (Int_t y = yMin; y <= yMax; y++) {
518       for (Int_t z = 0; z <= nBinsZ+1; z++) {
519         bin = x + nx * (y + ny * z);
520         Double_t val = histo3D->GetBinContent(bin);
521         projH->Fill(zAxis->GetBinCenter(z), val);
522         entries += val;
523       }
524     }
525   }
526   projH->SetEntries((Long64_t)entries);
527   return projH;
528 }
529
530 //_____________________________________________________________________
531 Long64_t AliTPCCalibKr::Merge(TCollection* list) {
532 // merge function 
533 //
534 cout << "Merge " << endl;
535
536 if (!list)
537 return 0;
538
539 if (list->IsEmpty())
540 return 1;
541
542 TIterator* iter = list->MakeIterator();
543 TObject* obj = 0;
544
545     iter->Reset();
546     Int_t count=0;
547         while((obj = iter->Next()) != 0)
548         {
549           AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
550           if (entry == 0) continue; 
551
552                 for(int i=0; i<72; ++i) { 
553                   if(IsCSide(i) == kTRUE && fCSide == kTRUE) { 
554                   ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
555                   } 
556
557                   if(IsCSide(i) == kFALSE && fASide == kTRUE) { 
558                     ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
559                   } 
560                 } 
561
562           count++;
563         }
564
565 return count;
566 }