07ab8473fa6fc3016953df7b46f053411ed247a0
[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 fills the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
52 // its data memebers.   
53 // 
54 // As the input it requires the tree with reconstructed Kr clusters (AliTPCclusterKr objects). 
55 // The AliTPCCalibKr objects containing an array of TH3F histograms are stored (by default) in the 
56 // ouptut (outHistFile.root) file.
57 //
58 // The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
59 // These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration 
60 // parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
61 // In addition the debugCalibKr.root file with debug information is created.
62 //
63 // Usage example:
64 //
65 // 1. Create outHistFile.root histogram file:
66 //
67 // -- Load libXrdClient.so if data on Xrd cluster e.g. (GSI)
68 // gSystem->Load("/usr/local/grid/XRootd/GSI/lib64/libXrdClient.so");
69 //
70 // -- Load toolkit
71 // gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
72 // gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
73 // AliXRDPROOFtoolkit tool;
74 //
75 // -- Make chain of files
76 // TChain * chain = tool.MakeChain("KrClusters.txt","Kr","",1000,0);
77 //
78 // -- Run AliTPCCalibKr task (Only TPC C side)
79 // AliTPCCalibKr *task = new AliTPCCalibKr;
80 // task->SetInputChain(chain);
81 // task->SetASide(kFALSE);
82 //
83 // task->Process();
84 // 
85 // 2. Analyse output histograms:
86 //
87 // TFile f("outHistFile.root");
88 // AliTPCCalibKr.Analyse();
89 //
90 // 3. See calibration parameters e.g.:
91 //
92 // TFile f("calibKr.root");
93 // spectrMean->GetCalROC(70)->GetValue(40,40);
94 // fitMean->GetCalROC(70)->GetValue(40,40);
95 //
96 // 4. See debug information e.g.:
97 //
98 // TFile f("debugCalibKr.root");
99 // .ls;
100 //
101 // -- Print calibKr TTree content 
102 // calibKr->Print();
103 //
104 // -- Draw calibKr TTree variables
105 // calibKr.Draw("fitMean");
106 //
107 //
108 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
109 //-----------------------------------------------------------------------------
110
111 ClassImp(AliTPCCalibKr)
112
113 AliTPCCalibKr::AliTPCCalibKr() : 
114   TObject(),
115   
116   bOutputHisto(kTRUE),
117   bASide(kTRUE),
118   bCSide(kTRUE),
119   fClusters(0),
120   fClustKr(0),
121   fTree(0),
122   fHistoKrArray(72)
123 {
124   //
125   // default constructor
126   //
127 }
128
129 //_____________________________________________________________________
130 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : 
131   TObject(pad),
132   
133   bOutputHisto(pad.bOutputHisto),
134   bASide(pad.bASide),
135   bCSide(pad.bCSide),
136   fClusters(pad.fClusters),
137   fClustKr(pad.fClustKr),
138   fTree(pad.fTree),
139   fHistoKrArray(72)
140 {
141   // copy constructor
142  
143   for (Int_t iSec = 0; iSec < 72; ++iSec) 
144   {
145     TH3F *hOld = pad.GetHistoKr(iSec);
146         if(hOld) {
147       TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); 
148       fHistoKrArray.AddAt(hNew,iSec);
149         }
150   }
151 }
152
153 //_____________________________________________________________________
154 AliTPCCalibKr::~AliTPCCalibKr() 
155 {
156   //
157   // destructor
158   //
159   if(fClustKr)  delete fClustKr; fClustKr = 0;
160   if(fClusters) delete fClusters; fClusters = 0;
161   if(fTree)     delete fTree; fTree = 0;
162   fHistoKrArray.Delete();
163 }
164
165 //_____________________________________________________________________
166 AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
167 {
168   // assignment operator
169
170   if (&source == this) return *this;
171   new (this) AliTPCCalibKr(source);
172
173   return *this;
174 }
175
176 //_____________________________________________________________________
177 void AliTPCCalibKr::Init()
178 {
179   // 
180   // init input tree and output histograms 
181   //
182
183   // set input tree
184   if(!fTree) { 
185    Printf("ERROR: Could not read chain from input");
186   }
187   else {
188    fTree->SetBranchStatus("*",1); 
189   }
190
191   // set branch address
192   fClusters = new TClonesArray("AliTPCclusterKr");
193
194   if(!fTree->GetBranch("fClusters")) {
195     Printf("ERROR: Could not get fClusters branch from input");
196   } else {
197    fTree->GetBranch("fClusters")->SetAddress(&fClusters);
198   }
199   
200   // create output TObjArray
201   fHistoKrArray.Clear();
202
203   // add histograms to the TObjArray
204   for(Int_t i=0; i<72; ++i) {
205     
206         // C - side
207         if( IsCSide(i) == kTRUE && bCSide == kTRUE) {
208       TH3F *hist = CreateHisto(i);
209       if(hist) fHistoKrArray.AddAt(hist,i);
210         }
211     
212         // A - side
213         if(IsCSide(i) == kFALSE && bASide == kTRUE) {
214       TH3F *hist = CreateHisto(i);
215       if(hist) fHistoKrArray.AddAt(hist,i);
216         }
217
218   }
219 }
220
221 //_____________________________________________________________________
222 Bool_t AliTPCCalibKr::ReadEntry(Int_t evt)
223 {
224   // 
225   // read entry from the tree
226   //
227   Long64_t centry = fTree->LoadTree(evt);
228   if(centry < 0) return kFALSE;
229
230   if(!fTree->GetBranch("fClusters")) 
231   {
232     Printf("ERROR: Could not get fClusters branch from input");
233         return kFALSE;
234   } else {
235    fTree->GetBranch("fClusters")->SetAddress(&fClusters);
236   }
237
238   fTree->GetEntry(evt);
239
240 return kTRUE;
241 }
242  
243 //_____________________________________________________________________
244 Bool_t AliTPCCalibKr::Process()
245 {
246   //
247   // process events 
248   // call event by event
249   //
250
251   // init tree
252   Init();
253
254   // get events
255   if(!fTree) return kFALSE;
256   Int_t nEvents = fTree->GetEntries();
257
258   // fill histograms 
259   for(Int_t i=0; i<nEvents; ++i)
260   {
261     if(ReadEntry(i) == kFALSE) return kFALSE;
262
263     if(!(i%10000)) cout << "evt: " << i << endl; 
264
265     // get TClonesArray entries
266     fClustKr = 0;
267     Int_t entries = fClusters->GetEntries();
268     for(Int_t j=0; j < entries; ++j)
269         {
270           fClustKr = (AliTPCclusterKr*)fClusters->At(j);
271
272       if(fClustKr) Update(fClustKr);
273           else return kFALSE;
274         }
275   }
276
277   // write output 
278   return Terminate();
279 }
280
281 //_____________________________________________________________________
282 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
283 {
284     //
285     // create new histogram
286         //
287     char name[256];
288         TH3F *h;
289
290     sprintf(name,"ADCcluster_ch%d",chamber);
291
292     if( IsIROC(chamber) == kTRUE ) 
293         {
294            h = new TH3F(name,name,63,0,63,100,0,100,150,100,3000);
295         } else {
296            h = new TH3F(name,name,96,0,96,100,0,100,150,100,3000);
297         }
298     h->SetXTitle("padrow");
299     h->SetYTitle("pad");
300     h->SetZTitle("fADC");
301
302 return h;
303 }
304
305 //_____________________________________________________________________
306 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
307 {
308 // check if IROCs
309 // returns kTRUE if IROCs and kFALSE if OROCs 
310
311   if(chamber>=0 && chamber<36) return kTRUE;
312
313 return kFALSE;
314 }
315
316 //_____________________________________________________________________
317 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
318 {
319 // check if C side
320 // returns kTRUE if C side and kFALSE if A side
321
322   if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
323
324 return kFALSE;
325 }
326  
327 //_____________________________________________________________________
328 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
329 {
330   //
331   // fill existing histograms
332   //
333
334   if (!Accept(cl)) return kFALSE;
335   TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
336   if(!h) return kFALSE;
337   
338   h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
339   
340   return kTRUE;
341 }
342
343 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
344   //
345   // cuts
346   //
347   /*
348     TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings - 
349     TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
350     TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
351     TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
352     TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
353     TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
354   */
355   //R0
356   if (cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
357   // R1
358   if (cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
359   //R2
360   if (cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
361   //R3
362   if (cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
363   //S1
364   if (cl->GetSize()>200) return kFALSE;
365   if (cl->GetSize()<6)  return kFALSE;
366   return kTRUE;
367
368 }
369
370
371
372 //_____________________________________________________________________
373 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
374 {
375   // get histograms from fHistoKrArray
376   return (TH3F*) fHistoKrArray.At(chamber);
377 }
378
379 //_____________________________________________________________________
380 Bool_t AliTPCCalibKr::Terminate() 
381 {
382   //
383   // store AliTPCCalibKr in the output file 
384   //
385   if(bOutputHisto) {
386     TFile *outFile = new TFile("outHistFile.root","RECREATE"); 
387    
388     if(outFile) 
389         {
390           outFile->cd();
391
392           for(int i=0; i<72; ++i) {
393              if( IsCSide(i) == kTRUE && bCSide == kTRUE)
394                printf("C side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
395
396              if( IsCSide(i) == kFALSE && bASide == kTRUE)
397                printf("A side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
398           }
399           this->Write();
400           outFile->Close();
401
402           return kTRUE;
403         }
404         else 
405           return kFALSE;
406   }
407
408 return kFALSE;
409 }
410  
411 //_____________________________________________________________________
412 void AliTPCCalibKr::Analyse() 
413 {
414   //
415   // analyse the histograms and extract krypton calibration parameters
416   //
417
418   // AliTPCCalPads that will contain the calibration parameters
419   AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
420   AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
421   AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
422   AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
423   AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
424   AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
425
426   // file stream for debugging purposes
427   TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
428
429   // if entries in spectrum less than minEntries, then the fit won't be performed
430   Int_t minEntries = 1; //300;
431
432   Double_t windowFrac = 0.12;
433   // the 3d histogram will be projected on the pads given by the following window size
434   // set the numbers to 0 if you want to do a pad-by-pad calibration
435   UInt_t rowRadius = 5;
436   UInt_t padRadius = 10;
437   // the step size by which pad and row are incremented is given by the following numbers
438   // set them to 1 if you want the finest granularity
439   UInt_t rowStep = 1;    // formerly: 2*rowRadius
440   UInt_t padStep = 1;    // formerly: 2*padRadius
441
442   for (Int_t chamber = 0; chamber < 72; chamber++) {
443     //if (chamber != 71) continue;
444     AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
445     
446     // Usually I would traverse each pad, take the spectrum for its neighbourhood and
447     // obtain the calibration parameters. This takes very long, so instead I assign the same
448     // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
449     UInt_t nRows = roc.GetNrows();
450     for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
451       UInt_t nPads = roc.GetNPads(iRow);
452       //if (iRow >= 10) break;
453       for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
454         //if (iPad >= 20) break;
455         TH3F* h = GetHistoKr(chamber);
456         if (!h) continue;
457         
458         // the 3d histogram will be projected on the pads given by the following bounds
459         // for rows and pads
460         Int_t rowLow = iRow - rowRadius;
461         UInt_t rowUp = iRow + rowRadius;
462         Int_t padLow = iPad - padRadius;
463         UInt_t padUp = iPad + padRadius;
464         // if window goes out of chamber
465         if (rowLow < 0) rowLow = 0;
466         if (rowUp >= nRows) rowUp = nRows - 1;
467         if (padLow < 0) padLow = 0;
468         if (padUp >= nPads) padUp = nPads - 1;
469
470         // project the histogram
471         //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
472         TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
473     
474         // get the number of entries in the spectrum
475         Double_t entries = projH->GetEntries();
476         if (entries < minEntries) { delete projH; continue; }
477         
478         // get the two calibration parameters mean of spectrum and RMS of spectrum
479         Double_t histMean = projH->GetMean();
480         Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
481     
482         // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
483         Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
484         Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
485         if (fitResult != 0) {
486           Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i.", chamber, rowLow, rowUp, padLow, padUp);
487           delete projH;
488           continue;
489         }
490     
491         // get the two calibration parameters mean of gauss fit and sigma of gauss fit
492         TF1* gausFit = projH->GetFunction("gaus");
493         Double_t fitMean = gausFit->GetParameter(1);
494         Double_t fitRMS = gausFit->GetParameter(2);
495         Int_t numberFitPoints = gausFit->GetNumberFitPoints();
496         if (numberFitPoints == 0) continue;
497         Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
498         delete gausFit;
499         delete projH;
500         if (fitMean <= 0) continue;
501         printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
502     
503         // write the calibration parameters for each pad that the 3d histogram was projected onto
504         // (with considering the step size) to the CalPads
505         // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
506         // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
507         for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
508           if (r < 0 || r >= (Int_t)nRows) continue;
509           UInt_t nPads = roc.GetNPads(r);
510           for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
511             if (p < 0 || p >= (Int_t)nPads) continue;
512             spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
513             spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
514             fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
515             fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
516             fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
517             entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
518
519             (*debugStream) << "calibKr" <<
520               "sector=" << chamber <<          // chamber number
521               "row=" << r <<                   // row number
522               "pad=" << p <<                   // pad number
523               "histMean=" << histMean <<       // mean of the spectrum
524               "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
525               "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
526               "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
527               "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
528               "entries=" << entries <<         // number of entries for the spectrum
529               "\n";
530           }
531         }
532       }
533     }
534   }
535
536   TFile f("calibKr.root", "recreate");
537   spectrMeanCalPad->Write();
538   spectrRMSCalPad->Write();
539   fitMeanCalPad->Write();
540   fitRMSCalPad->Write();
541   fitNormChi2CalPad->Write();
542   entriesCalPad->Write();
543   f.Close();
544   delete spectrMeanCalPad;
545   delete spectrRMSCalPad;
546   delete fitMeanCalPad;
547   delete fitRMSCalPad;
548   delete fitNormChi2CalPad;
549   delete entriesCalPad;
550   delete debugStream;
551 }
552
553 //_____________________________________________________________________
554 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
555 {
556   // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
557   // replaces TH3F::ProjectZ() to gain more speed
558
559   TAxis* xAxis = histo3D->GetXaxis();
560   TAxis* yAxis = histo3D->GetYaxis();
561   TAxis* zAxis = histo3D->GetZaxis();
562   Double_t zMinVal = zAxis->GetXmin();
563   Double_t zMaxVal = zAxis->GetXmax();
564   
565   Int_t nBinsZ = zAxis->GetNbins();
566   TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
567
568   Int_t nx = xAxis->GetNbins()+2;
569   Int_t ny = yAxis->GetNbins()+2;
570   Int_t bin = 0;
571   Double_t entries = 0.;
572   for (Int_t x = xMin; x <= xMax; x++) {
573     for (Int_t y = yMin; y <= yMax; y++) {
574       for (Int_t z = 0; z <= nBinsZ+1; z++) {
575         bin = x + nx * (y + ny * z);
576         Double_t val = histo3D->GetBinContent(bin);
577         projH->Fill(zAxis->GetBinCenter(z), val);
578         entries += val;
579       }
580     }
581   }
582   projH->SetEntries((Long64_t)entries);
583   return projH;
584 }