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 fills the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
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.
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.
65 // 1. Create outHistFile.root histogram file:
67 // -- Load libXrdClient.so if data on Xrd cluster e.g. (GSI)
68 // gSystem->Load("/usr/local/grid/XRootd/GSI/lib64/libXrdClient.so");
71 // gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
72 // gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
73 // AliXRDPROOFtoolkit tool;
75 // -- Make chain of files
76 // TChain * chain = tool.MakeChain("KrClusters.txt","Kr","",1000,0);
78 // -- Run AliTPCCalibKr task (Only TPC C side)
79 // AliTPCCalibKr *task = new AliTPCCalibKr;
80 // task->SetInputChain(chain);
81 // task->SetASide(kFALSE);
85 // 2. Analyse output histograms:
87 // TFile f("outHistFile.root");
88 // AliTPCCalibKr.Analyse();
90 // 3. See calibration parameters e.g.:
92 // TFile f("calibKr.root");
93 // spectrMean->GetCalROC(70)->GetValue(40,40);
94 // fitMean->GetCalROC(70)->GetValue(40,40);
96 // 4. See debug information e.g.:
98 // TFile f("debugCalibKr.root");
101 // -- Print calibKr TTree content
104 // -- Draw calibKr TTree variables
105 // calibKr.Draw("fitMean");
108 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
109 //-----------------------------------------------------------------------------
111 ClassImp(AliTPCCalibKr)
113 AliTPCCalibKr::AliTPCCalibKr() :
125 // default constructor
129 //_____________________________________________________________________
130 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
133 bOutputHisto(pad.bOutputHisto),
136 fClusters(pad.fClusters),
137 fClustKr(pad.fClustKr),
143 for (Int_t iSec = 0; iSec < 72; ++iSec)
145 TH3F *hOld = pad.GetHistoKr(iSec);
147 TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
148 fHistoKrArray.AddAt(hNew,iSec);
153 //_____________________________________________________________________
154 AliTPCCalibKr::~AliTPCCalibKr()
159 if(fClustKr) delete fClustKr; fClustKr = 0;
160 if(fClusters) delete fClusters; fClusters = 0;
161 if(fTree) delete fTree; fTree = 0;
162 fHistoKrArray.Delete();
165 //_____________________________________________________________________
166 AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
168 // assignment operator
170 if (&source == this) return *this;
171 new (this) AliTPCCalibKr(source);
176 //_____________________________________________________________________
177 void AliTPCCalibKr::Init()
180 // init input tree and output histograms
185 Printf("ERROR: Could not read chain from input");
188 fTree->SetBranchStatus("*",1);
191 // set branch address
192 fClusters = new TClonesArray("AliTPCclusterKr");
194 if(!fTree->GetBranch("fClusters")) {
195 Printf("ERROR: Could not get fClusters branch from input");
197 fTree->GetBranch("fClusters")->SetAddress(&fClusters);
200 // create output TObjArray
201 fHistoKrArray.Clear();
203 // add histograms to the TObjArray
204 for(Int_t i=0; i<72; ++i) {
207 if( IsCSide(i) == kTRUE && bCSide == kTRUE) {
208 TH3F *hist = CreateHisto(i);
209 if(hist) fHistoKrArray.AddAt(hist,i);
213 if(IsCSide(i) == kFALSE && bASide == kTRUE) {
214 TH3F *hist = CreateHisto(i);
215 if(hist) fHistoKrArray.AddAt(hist,i);
221 //_____________________________________________________________________
222 Bool_t AliTPCCalibKr::ReadEntry(Int_t evt)
225 // read entry from the tree
227 Long64_t centry = fTree->LoadTree(evt);
228 if(centry < 0) return kFALSE;
230 if(!fTree->GetBranch("fClusters"))
232 Printf("ERROR: Could not get fClusters branch from input");
235 fTree->GetBranch("fClusters")->SetAddress(&fClusters);
238 fTree->GetEntry(evt);
243 //_____________________________________________________________________
244 Bool_t AliTPCCalibKr::Process()
248 // call event by event
255 if(!fTree) return kFALSE;
256 Int_t nEvents = fTree->GetEntries();
259 for(Int_t i=0; i<nEvents; ++i)
261 if(ReadEntry(i) == kFALSE) return kFALSE;
263 if(!(i%10000)) cout << "evt: " << i << endl;
265 // get TClonesArray entries
267 Int_t entries = fClusters->GetEntries();
268 for(Int_t j=0; j < entries; ++j)
270 fClustKr = (AliTPCclusterKr*)fClusters->At(j);
272 if(fClustKr) Update(fClustKr);
281 //_____________________________________________________________________
282 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
285 // create new histogram
290 sprintf(name,"ADCcluster_ch%d",chamber);
292 if( IsIROC(chamber) == kTRUE )
294 h = new TH3F(name,name,63,0,63,100,0,100,150,100,3000);
296 h = new TH3F(name,name,96,0,96,100,0,100,150,100,3000);
298 h->SetXTitle("padrow");
300 h->SetZTitle("fADC");
305 //_____________________________________________________________________
306 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
309 // returns kTRUE if IROCs and kFALSE if OROCs
311 if(chamber>=0 && chamber<36) return kTRUE;
316 //_____________________________________________________________________
317 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
320 // returns kTRUE if C side and kFALSE if A side
322 if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
327 //_____________________________________________________________________
328 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
331 // fill existing histograms
333 TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
334 if(!h) return kFALSE;
337 h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
342 //_____________________________________________________________________
343 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
345 // get histograms from fHistoKrArray
346 return (TH3F*) fHistoKrArray.At(chamber);
349 //_____________________________________________________________________
350 Bool_t AliTPCCalibKr::Terminate()
353 // store AliTPCCalibKr in the output file
356 TFile *outFile = new TFile("outHistFile.root","RECREATE");
362 for(int i=0; i<72; ++i) {
363 if( IsCSide(i) == kTRUE && bCSide == kTRUE)
364 printf("C side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
366 if( IsCSide(i) == kFALSE && bASide == kTRUE)
367 printf("A side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
381 //_____________________________________________________________________
382 void AliTPCCalibKr::Analyse()
385 // analyse the histograms and extract krypton calibration parameters
388 // AliTPCCalPads that will contain the calibration parameters
389 AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
390 AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
391 AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
392 AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
393 AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
394 AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
396 // file stream for debugging purposes
397 TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
399 // if entries in spectrum less than minEntries, then the fit won't be performed
400 Int_t minEntries = 1; //300;
402 Double_t windowFrac = 0.12;
403 // the 3d histogram will be projected on the pads given by the following window size
404 // set the numbers to 0 if you want to do a pad-by-pad calibration
405 UInt_t rowRadius = 5;
406 UInt_t padRadius = 10;
407 // the step size by which pad and row are incremented is given by the following numbers
408 // set them to 1 if you want the finest granularity
409 UInt_t rowStep = 1; // formerly: 2*rowRadius
410 UInt_t padStep = 1; // formerly: 2*padRadius
412 for (Int_t chamber = 0; chamber < 72; chamber++) {
413 //if (chamber != 71) continue;
414 AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()
416 // Usually I would traverse each pad, take the spectrum for its neighbourhood and
417 // obtain the calibration parameters. This takes very long, so instead I assign the same
418 // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
419 UInt_t nRows = roc.GetNrows();
420 for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
421 UInt_t nPads = roc.GetNPads(iRow);
422 //if (iRow >= 10) break;
423 for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
424 //if (iPad >= 20) break;
425 TH3F* h = GetHistoKr(chamber);
428 // the 3d histogram will be projected on the pads given by the following bounds
430 Int_t rowLow = iRow - rowRadius;
431 UInt_t rowUp = iRow + rowRadius;
432 Int_t padLow = iPad - padRadius;
433 UInt_t padUp = iPad + padRadius;
434 // if window goes out of chamber
435 if (rowLow < 0) rowLow = 0;
436 if (rowUp >= nRows) rowUp = nRows - 1;
437 if (padLow < 0) padLow = 0;
438 if (padUp >= nPads) padUp = nPads - 1;
440 // project the histogram
441 //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
442 TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
444 // get the number of entries in the spectrum
445 Double_t entries = projH->GetEntries();
446 if (entries < minEntries) { delete projH; continue; }
448 // get the two calibration parameters mean of spectrum and RMS of spectrum
449 Double_t histMean = projH->GetMean();
450 Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
452 // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
453 Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
454 Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
455 if (fitResult != 0) {
456 Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i.", chamber, rowLow, rowUp, padLow, padUp);
461 // get the two calibration parameters mean of gauss fit and sigma of gauss fit
462 TF1* gausFit = projH->GetFunction("gaus");
463 Double_t fitMean = gausFit->GetParameter(1);
464 Double_t fitRMS = gausFit->GetParameter(2);
465 Int_t numberFitPoints = gausFit->GetNumberFitPoints();
466 if (numberFitPoints == 0) continue;
467 Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
470 if (fitMean <= 0) continue;
471 printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
473 // write the calibration parameters for each pad that the 3d histogram was projected onto
474 // (with considering the step size) to the CalPads
475 // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
476 // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
477 for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
478 if (r < 0 || r >= (Int_t)nRows) continue;
479 UInt_t nPads = roc.GetNPads(r);
480 for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
481 if (p < 0 || p >= (Int_t)nPads) continue;
482 spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
483 spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
484 fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
485 fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
486 fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
487 entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
489 (*debugStream) << "calibKr" <<
490 "sector=" << chamber << // chamber number
491 "row=" << r << // row number
492 "pad=" << p << // pad number
493 "histMean=" << histMean << // mean of the spectrum
494 "histRMS=" << histRMS << // RMS of the spectrum divided by the mean
495 "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak
496 "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak
497 "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit
498 "entries=" << entries << // number of entries for the spectrum
506 TFile f("calibKr.root", "recreate");
507 spectrMeanCalPad->Write();
508 spectrRMSCalPad->Write();
509 fitMeanCalPad->Write();
510 fitRMSCalPad->Write();
511 fitNormChi2CalPad->Write();
512 entriesCalPad->Write();
514 delete spectrMeanCalPad;
515 delete spectrRMSCalPad;
516 delete fitMeanCalPad;
518 delete fitNormChi2CalPad;
519 delete entriesCalPad;
523 //_____________________________________________________________________
524 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
526 // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
527 // replaces TH3F::ProjectZ() to gain more speed
529 TAxis* xAxis = histo3D->GetXaxis();
530 TAxis* yAxis = histo3D->GetYaxis();
531 TAxis* zAxis = histo3D->GetZaxis();
532 Double_t zMinVal = zAxis->GetXmin();
533 Double_t zMaxVal = zAxis->GetXmax();
535 Int_t nBinsZ = zAxis->GetNbins();
536 TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
538 Int_t nx = xAxis->GetNbins()+2;
539 Int_t ny = yAxis->GetNbins()+2;
541 Double_t entries = 0.;
542 for (Int_t x = xMin; x <= xMax; x++) {
543 for (Int_t y = yMin; y <= yMax; y++) {
544 for (Int_t z = 0; z <= nBinsZ+1; z++) {
545 bin = x + nx * (y + ny * z);
546 Double_t val = histo3D->GetBinContent(bin);
547 projH->Fill(zAxis->GetBinCenter(z), val);
552 projH->SetEntries((Long64_t)entries);