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() :
101 // default constructor
104 // TObjArray with histograms
105 fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
106 fHistoKrArray.Clear();
112 //_____________________________________________________________________
113 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
122 // TObjArray with histograms
123 fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
124 fHistoKrArray.Clear();
126 for (Int_t iSec = 0; iSec < 72; ++iSec)
128 TH3F *hOld = pad.GetHistoKr(iSec);
130 TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
131 fHistoKrArray.AddAt(hNew,iSec);
136 //_____________________________________________________________________
137 AliTPCCalibKr::~AliTPCCalibKr()
143 // for (Int_t iSec = 0; iSec < 72; ++iSec) {
144 // if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
146 fHistoKrArray.Delete();
149 //_____________________________________________________________________
150 AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
152 // assignment operator
154 if (&source == this) return *this;
155 new (this) AliTPCCalibKr(source);
160 //_____________________________________________________________________
161 void AliTPCCalibKr::Init()
164 // init output histograms
167 // add histograms to the TObjArray
168 for(Int_t i=0; i<72; ++i) {
171 if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
172 TH3F *hist = CreateHisto(i);
173 if(hist) fHistoKrArray.AddAt(hist,i);
177 if(IsCSide(i) == kFALSE && fASide == kTRUE) {
178 TH3F *hist = CreateHisto(i);
179 if(hist) fHistoKrArray.AddAt(hist,i);
184 //_____________________________________________________________________
185 Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
189 // call event by event
192 if(cluster) Update(cluster);
198 //_____________________________________________________________________
199 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
202 // create new histogram
207 sprintf(name,"ADCcluster_ch%d",chamber);
209 if( IsIROC(chamber) == kTRUE )
211 h = new TH3F(name,name,63,0,63,108,0,108,200,100,2500);
213 h = new TH3F(name,name,96,0,96,140,0,140,200,100,1700);
215 h->SetXTitle("padrow");
217 h->SetZTitle("fADC");
222 //_____________________________________________________________________
223 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
226 // returns kTRUE if IROCs and kFALSE if OROCs
228 if(chamber>=0 && chamber<36) return kTRUE;
233 //_____________________________________________________________________
234 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
237 // returns kTRUE if C side and kFALSE if A side
239 if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
244 //_____________________________________________________________________
245 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
248 // fill existing histograms
251 if (!Accept(cl)) return kFALSE;
252 TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
253 if(!h) return kFALSE;
255 h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
260 //_____________________________________________________________________
261 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){
266 TCut cutR0("cutR0","fADCcluster/fSize<200"); // adjust it according v seetings -
267 TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
268 TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4"); // digital noise removal
269 TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
270 TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
271 TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
274 if ((float)cl->GetADCcluster()/ cl->GetSize() >200) return kFALSE;
276 if ((float)cl->GetADCcluster()/ cl->GetSize() <7) return kFALSE;
278 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4) return kFALSE;
280 if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
282 if (cl->GetSize()>200) return kFALSE;
283 if (cl->GetSize()<6) return kFALSE;
288 //_____________________________________________________________________
289 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
291 // get histograms from fHistoKrArray
292 return (TH3F*) fHistoKrArray.At(chamber);
295 //_____________________________________________________________________
296 void AliTPCCalibKr::Analyse()
299 // analyse the histograms and extract krypton calibration parameters
302 // AliTPCCalPads that will contain the calibration parameters
303 AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
304 AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
305 AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
306 AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
307 AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
308 AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
310 // file stream for debugging purposes
311 TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
313 // if entries in spectrum less than minEntries, then the fit won't be performed
314 Int_t minEntries = 1; //300;
316 Double_t windowFrac = 0.12;
317 // the 3d histogram will be projected on the pads given by the following window size
318 // set the numbers to 0 if you want to do a pad-by-pad calibration
319 UInt_t rowRadius = 2;
320 UInt_t padRadius = 4;
322 // the step size by which pad and row are incremented is given by the following numbers
323 // set them to 1 if you want the finest granularity
324 UInt_t rowStep = 1; // formerly: 2*rowRadius
325 UInt_t padStep = 1; // formerly: 2*padRadius
327 for (Int_t chamber = 0; chamber < 72; chamber++) {
328 //if (chamber != 71) continue;
329 AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()
331 // Usually I would traverse each pad, take the spectrum for its neighbourhood and
332 // obtain the calibration parameters. This takes very long, so instead I assign the same
333 // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
334 UInt_t nRows = roc.GetNrows();
335 for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
336 UInt_t nPads = roc.GetNPads(iRow);
337 //if (iRow >= 10) break;
338 for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
339 //if (iPad >= 20) break;
340 TH3F* h = GetHistoKr(chamber);
343 // the 3d histogram will be projected on the pads given by the following bounds
345 Int_t rowLow = iRow - rowRadius;
346 UInt_t rowUp = iRow + rowRadius;
347 Int_t padLow = iPad - padRadius;
348 UInt_t padUp = iPad + padRadius;
349 // if window goes out of chamber
350 if (rowLow < 0) rowLow = 0;
351 if (rowUp >= nRows) rowUp = nRows - 1;
352 if (padLow < 0) padLow = 0;
353 if (padUp >= nPads) padUp = nPads - 1;
355 // project the histogram
356 //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
357 TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
359 // get the number of entries in the spectrum
360 Double_t entries = projH->GetEntries();
361 if (entries < minEntries) { delete projH; continue; }
363 // get the two calibration parameters mean of spectrum and RMS of spectrum
364 Double_t histMean = projH->GetMean();
365 Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
367 // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
368 Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
369 Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
370 Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
371 Double_t integCharge = projH->Integral(minBin,maxBin);
373 Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
375 if (fitResult != 0) {
376 Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
377 //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, %f\n", chamber, iRow, iPad, entries, maxEntries);
383 // get the two calibration parameters mean of gauss fit and sigma of gauss fit
384 TF1* gausFit = projH->GetFunction("gaus");
385 Double_t fitMean = gausFit->GetParameter(1);
386 Double_t fitRMS = gausFit->GetParameter(2);
387 Int_t numberFitPoints = gausFit->GetNumberFitPoints();
388 if (numberFitPoints == 0) continue;
389 Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
392 if (fitMean <= 0) continue;
393 //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
395 // write the calibration parameters for each pad that the 3d histogram was projected onto
396 // (with considering the step size) to the CalPads
397 // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
398 // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
399 for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
400 if (r < 0 || r >= (Int_t)nRows) continue;
401 UInt_t nPads = roc.GetNPads(r);
402 for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
403 if (p < 0 || p >= (Int_t)nPads) continue;
404 spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
405 spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
406 fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
407 fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
408 fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
409 entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
411 (*debugStream) << "calibKr" <<
412 "sector=" << chamber << // chamber number
413 "row=" << r << // row number
414 "pad=" << p << // pad number
415 "histMean=" << histMean << // mean of the spectrum
416 "histRMS=" << histRMS << // RMS of the spectrum divided by the mean
417 "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak
418 "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak
419 "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit
420 "entries=" << entries << // number of entries for the spectrum
428 TFile f("calibKr.root", "recreate");
429 spectrMeanCalPad->Write();
430 spectrRMSCalPad->Write();
431 fitMeanCalPad->Write();
432 fitRMSCalPad->Write();
433 fitNormChi2CalPad->Write();
434 entriesCalPad->Write();
436 delete spectrMeanCalPad;
437 delete spectrRMSCalPad;
438 delete fitMeanCalPad;
440 delete fitNormChi2CalPad;
441 delete entriesCalPad;
445 //_____________________________________________________________________
446 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
448 // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
449 // replaces TH3F::ProjectZ() to gain more speed
451 TAxis* xAxis = histo3D->GetXaxis();
452 TAxis* yAxis = histo3D->GetYaxis();
453 TAxis* zAxis = histo3D->GetZaxis();
454 Double_t zMinVal = zAxis->GetXmin();
455 Double_t zMaxVal = zAxis->GetXmax();
457 Int_t nBinsZ = zAxis->GetNbins();
458 TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
460 Int_t nx = xAxis->GetNbins()+2;
461 Int_t ny = yAxis->GetNbins()+2;
463 Double_t entries = 0.;
464 for (Int_t x = xMin; x <= xMax; x++) {
465 for (Int_t y = yMin; y <= yMax; y++) {
466 for (Int_t z = 0; z <= nBinsZ+1; z++) {
467 bin = x + nx * (y + ny * z);
468 Double_t val = histo3D->GetBinContent(bin);
469 projH->Fill(zAxis->GetBinCenter(z), val);
474 projH->SetEntries((Long64_t)entries);
478 //_____________________________________________________________________
479 Long64_t AliTPCCalibKr::Merge(TCollection* list) {
482 cout << "Merge " << endl;
490 TIterator* iter = list->MakeIterator();
495 while((obj = iter->Next()) != 0)
497 AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
498 if (entry == 0) continue;
500 for(int i=0; i<72; ++i) {
501 if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
502 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
505 if(IsCSide(i) == kFALSE && fASide == kTRUE) {
506 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );