1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
26 #include <TDirectory.h>
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"
45 #include "AliTPCCalibKr.h"
47 //----------------------------------------------------------------------------
48 // The AliTPCCalibKr class description (TPC Kr calibration).
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()).
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.
65 // 1. Analyse output histograms:
66 TFile f("outHistFile.root");
67 AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr")
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);
75 // 3. See debug information e.g.:
76 TFile f("debugCalibKr.root");
79 // -- Print calibKr TTree content
82 // -- Draw calibKr TTree variables
83 calibKr.Draw("fitMean");
89 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
90 //-----------------------------------------------------------------------------
92 ClassImp(AliTPCCalibKr)
94 AliTPCCalibKr::AliTPCCalibKr() :
99 fADCOverClustSizeMin(0.0),
100 fADCOverClustSizeMax(1.0e9),
101 fMaxADCOverClustADCMin(0.0),
102 fMaxADCOverClustADCMax(1.0e9),
109 // default constructor
112 // TObjArray with histograms
113 fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
114 fHistoKrArray.Clear();
117 SetADCOverClustSizeRange(7,200);
118 SetMaxADCOverClustADCRange(0.01,0.4);
119 SetTimeRange(200,600);
120 SetClustSizeRange(6,200);
126 //_____________________________________________________________________
127 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
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)
144 // TObjArray with histograms
145 fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
146 fHistoKrArray.Clear();
148 for (Int_t iSec = 0; iSec < 72; ++iSec)
150 TH3F *hOld = pad.GetHistoKr(iSec);
152 TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
153 fHistoKrArray.AddAt(hNew,iSec);
158 //_____________________________________________________________________
159 AliTPCCalibKr::~AliTPCCalibKr()
165 // for (Int_t iSec = 0; iSec < 72; ++iSec) {
166 // if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
168 fHistoKrArray.Delete();
171 //_____________________________________________________________________
172 AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
174 // assignment operator
176 if (&source == this) return *this;
177 new (this) AliTPCCalibKr(source);
182 //_____________________________________________________________________
183 void AliTPCCalibKr::Init()
186 // init output histograms
189 // add histograms to the TObjArray
190 for(Int_t i=0; i<72; ++i) {
193 if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
194 TH3F *hist = CreateHisto(i);
195 if(hist) fHistoKrArray.AddAt(hist,i);
199 if(IsCSide(i) == kFALSE && fASide == kTRUE) {
200 TH3F *hist = CreateHisto(i);
201 if(hist) fHistoKrArray.AddAt(hist,i);
206 //_____________________________________________________________________
207 Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
211 // call event by event
214 if(cluster) Update(cluster);
220 //_____________________________________________________________________
221 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
224 // create new histogram
229 sprintf(name,"ADCcluster_ch%d",chamber);
231 if( IsIROC(chamber) == kTRUE )
233 h = new TH3F(name,name,63,0,63,108,0,108,200,100,2500);
235 h = new TH3F(name,name,96,0,96,140,0,140,200,100,1700);
237 h->SetXTitle("padrow");
239 h->SetZTitle("fADC");
244 //_____________________________________________________________________
245 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
248 // returns kTRUE if IROCs and kFALSE if OROCs
250 if(chamber>=0 && chamber<36) return kTRUE;
255 //_____________________________________________________________________
256 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
259 // returns kTRUE if C side and kFALSE if A side
261 if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
266 //_____________________________________________________________________
267 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
270 // fill existing histograms
273 if (!Accept(cl)) return kFALSE;
274 TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
275 if(!h) return kFALSE;
277 h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
282 //_____________________________________________________________________
283 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){
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;
299 if ((float)cl->GetADCcluster()/ cl->GetSize() >200) return kFALSE;
301 if ((float)cl->GetADCcluster()/ cl->GetSize() <7) return kFALSE;
303 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4) return kFALSE;
305 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
307 if (cl->GetMax().GetTime() < 200) return kFALSE;
309 if (cl->GetMax().GetTime() > 600) return kFALSE;
311 if (cl->GetSize()>200) return kFALSE;
312 if (cl->GetSize()<6) return kFALSE;
314 SetADCOverClustSizeRange(7,200);
315 SetMaxADCOverClustADCRange(0.01,0.4);
316 SetTimeRange(200,600);
317 SetClustSizeRange(6,200);
321 if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax) return kFALSE;
323 if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin) return kFALSE;
325 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax) return kFALSE;
327 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
329 if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
331 if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
333 if (cl->GetSize()>fClustSizeMax) return kFALSE;
334 if (cl->GetSize()<fClustSizeMin) return kFALSE;
340 //_____________________________________________________________________
341 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
343 // get histograms from fHistoKrArray
344 return (TH3F*) fHistoKrArray.At(chamber);
347 //_____________________________________________________________________
348 void AliTPCCalibKr::Analyse()
351 // analyse the histograms and extract krypton calibration parameters
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");
362 // file stream for debugging purposes
363 TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
365 // if entries in spectrum less than minEntries, then the fit won't be performed
366 Int_t minEntries = 1; //300;
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;
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
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()
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);
395 // the 3d histogram will be projected on the pads given by the following bounds
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;
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);
411 // get the number of entries in the spectrum
412 Double_t entries = projH->GetEntries();
413 if (entries < minEntries) { delete projH; continue; }
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.;
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);
425 Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
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);
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;
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);
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 nPadsR = 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)nPadsR) 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);
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
480 TFile f("calibKr.root", "recreate");
481 spectrMeanCalPad->Write();
482 spectrRMSCalPad->Write();
483 fitMeanCalPad->Write();
484 fitRMSCalPad->Write();
485 fitNormChi2CalPad->Write();
486 entriesCalPad->Write();
488 delete spectrMeanCalPad;
489 delete spectrRMSCalPad;
490 delete fitMeanCalPad;
492 delete fitNormChi2CalPad;
493 delete entriesCalPad;
497 //_____________________________________________________________________
498 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
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
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();
509 Int_t nBinsZ = zAxis->GetNbins();
510 TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
512 Int_t nx = xAxis->GetNbins()+2;
513 Int_t ny = yAxis->GetNbins()+2;
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);
526 projH->SetEntries((Long64_t)entries);
530 //_____________________________________________________________________
531 Long64_t AliTPCCalibKr::Merge(TCollection* list) {
534 cout << "Merge " << endl;
542 TIterator* iter = list->MakeIterator();
547 while((obj = iter->Next()) != 0)
549 AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
550 if (entry == 0) continue;
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)) );
557 if(IsCSide(i) == kFALSE && fASide == kTRUE) {
558 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );