-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-
-//Root includes
-#include <TH1F.h>
-#include <TH1D.h>
-#include <TH2F.h>
-#include <TH3F.h>
-#include <TString.h>
-#include <TMath.h>
-#include <TF1.h>
-#include <TRandom.h>
-#include <TDirectory.h>
-#include <TFile.h>
-#include <TAxis.h>
-//AliRoot includes
-#include "AliRawReader.h"
-#include "AliRawReaderRoot.h"
-#include "AliRawReaderDate.h"
-#include "AliTPCRawStream.h"
-#include "AliTPCCalROC.h"
-#include "AliTPCCalPad.h"
-#include "AliTPCROC.h"
-#include "AliMathBase.h"
-#include "TTreeStream.h"
-#include "AliTPCRawStreamFast.h"
-
-//date
-#include "event.h"
-
-//header file
-#include "AliTPCCalibKr.h"
-
-//----------------------------------------------------------------------------
-// The AliTPCCalibKr class description (TPC Kr calibration).
-//
-//
-// The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
-// its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).
-//
-// The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
-// These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration
-// parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
-// In addition the debugCalibKr.root file with debug information is created.
-//
-
-/*
-
-// Usage example:
-//
-
-// 1. Analyse output histograms:
-TFile f("outHistFile.root");
-AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr")
-obj->Analyse();
-
-// 2. See calibration parameters e.g.:
-TFile f("calibKr.root");
-spectrMean->GetCalROC(70)->GetValue(40,40);
-fitMean->GetCalROC(70)->GetValue(40,40);
-
-// 3. See debug information e.g.:
-TFile f("debugCalibKr.root");
-.ls;
-
-// -- Print calibKr TTree content
-calibKr->Print();
-
-// -- Draw calibKr TTree variables
-calibKr.Draw("fitMean");
-
-*/
-
-
-//
-// Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
-//-----------------------------------------------------------------------------
-
-ClassImp(AliTPCCalibKr)
-
-AliTPCCalibKr::AliTPCCalibKr() :
- TObject(),
-
- bASide(kTRUE),
- bCSide(kTRUE),
- fHistoKrArray(72)
-{
- //
- // default constructor
- //
-
- // init histograms
- Init();
-}
-
-//_____________________________________________________________________
-AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
- TObject(pad),
-
- bASide(pad.bASide),
- bCSide(pad.bCSide),
- fHistoKrArray(72)
-{
- // copy constructor
-
- for (Int_t iSec = 0; iSec < 72; ++iSec)
- {
- TH3F *hOld = pad.GetHistoKr(iSec);
- if(hOld) {
- TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
- fHistoKrArray.AddAt(hNew,iSec);
- }
- }
-}
-
-//_____________________________________________________________________
-AliTPCCalibKr::~AliTPCCalibKr()
-{
- //
- // destructor
- //
-
- for (Int_t iSec = 0; iSec < 72; ++iSec)
- {
- TH3F *hNew = (TH3F*)fHistoKrArray.At(iSec);
- if(hNew) delete hNew; hNew = 0;
- }
- fHistoKrArray.Delete();
-}
-
-//_____________________________________________________________________
-AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
-{
- // assignment operator
-
- if (&source == this) return *this;
- new (this) AliTPCCalibKr(source);
-
- return *this;
-}
-
-//_____________________________________________________________________
-void AliTPCCalibKr::Init()
-{
- //
- // init input tree and output histograms
- //
-
- // create output TObjArray
- fHistoKrArray.Clear();
-
- // add histograms to the TObjArray
- for(Int_t i=0; i<72; ++i) {
-
- // C - side
- if( IsCSide(i) == kTRUE && bCSide == kTRUE) {
- TH3F *hist = CreateHisto(i);
- if(hist) fHistoKrArray.AddAt(hist,i);
- }
-
- // A - side
- if(IsCSide(i) == kFALSE && bASide == kTRUE) {
- TH3F *hist = CreateHisto(i);
- if(hist) fHistoKrArray.AddAt(hist,i);
- }
- }
-}
-
-//_____________________________________________________________________
-Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
-{
- //
- // process events
- // call event by event
- //
-
- if(cluster) Update(cluster);
- else return kFALSE;
-}
-
-//_____________________________________________________________________
-TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
-{
- //
- // create new histogram
- //
- char name[256];
- TH3F *h;
-
- sprintf(name,"ADCcluster_ch%d",chamber);
-
- if( IsIROC(chamber) == kTRUE )
- {
- h = new TH3F(name,name,63,0,63,100,0,100,150,100,3000);
- } else {
- h = new TH3F(name,name,96,0,96,100,0,100,150,100,3000);
- }
- h->SetXTitle("padrow");
- h->SetYTitle("pad");
- h->SetZTitle("fADC");
-
-return h;
-}
-
-//_____________________________________________________________________
-Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
-{
-// check if IROCs
-// returns kTRUE if IROCs and kFALSE if OROCs
-
- if(chamber>=0 && chamber<36) return kTRUE;
-
-return kFALSE;
-}
-
-//_____________________________________________________________________
-Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
-{
-// check if C side
-// returns kTRUE if C side and kFALSE if A side
-
- if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
-
-return kFALSE;
-}
-
-//_____________________________________________________________________
-Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
-{
- //
- // fill existing histograms
- //
-
- if (!Accept(cl)) return kFALSE;
- TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
- if(!h) return kFALSE;
-
- h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
-
- return kTRUE;
-}
-
-//_____________________________________________________________________
-Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){
- //
- // cuts
- //
- /*
- TCut cutR0("cutR0","fADCcluster/fSize<200"); // adjust it according v seetings -
- TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
- TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4"); // digital noise removal
- TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
- TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
- TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
- */
- //R0
- if (cl->GetADCcluster()/ Float_t(cl->GetSize()) >200) return kFALSE;
- // R1
- if (cl->GetADCcluster()/ Float_t(cl->GetSize()) <7) return kFALSE;
- //R2
- if (cl->GetMax().GetAdc()/ Float_t(cl->GetADCcluster()) >0.4) return kFALSE;
- //R3
- if (cl->GetMax().GetAdc()/ Float_t(cl->GetADCcluster()) <0.01) return kFALSE;
- //S1
- if (cl->GetSize()>200) return kFALSE;
- if (cl->GetSize()<6) return kFALSE;
- return kTRUE;
-
-}
-
-//_____________________________________________________________________
-TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
-{
- // get histograms from fHistoKrArray
- return (TH3F*) fHistoKrArray.At(chamber);
-}
-
-//_____________________________________________________________________
-void AliTPCCalibKr::Analyse()
-{
- //
- // analyse the histograms and extract krypton calibration parameters
- //
-
- // AliTPCCalPads that will contain the calibration parameters
- AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
- AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
- AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
- AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
- AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
- AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
-
- // file stream for debugging purposes
- TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
-
- // if entries in spectrum less than minEntries, then the fit won't be performed
- Int_t minEntries = 1; //300;
-
- Double_t windowFrac = 0.12;
- // the 3d histogram will be projected on the pads given by the following window size
- // set the numbers to 0 if you want to do a pad-by-pad calibration
- UInt_t rowRadius = 5;
- UInt_t padRadius = 10;
- // the step size by which pad and row are incremented is given by the following numbers
- // set them to 1 if you want the finest granularity
- UInt_t rowStep = 1; // formerly: 2*rowRadius
- UInt_t padStep = 1; // formerly: 2*padRadius
-
- for (Int_t chamber = 0; chamber < 72; chamber++) {
- //if (chamber != 71) continue;
- AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()
-
- // Usually I would traverse each pad, take the spectrum for its neighbourhood and
- // obtain the calibration parameters. This takes very long, so instead I assign the same
- // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
- UInt_t nRows = roc.GetNrows();
- for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
- UInt_t nPads = roc.GetNPads(iRow);
- //if (iRow >= 10) break;
- for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
- //if (iPad >= 20) break;
- TH3F* h = GetHistoKr(chamber);
- if (!h) continue;
-
- // the 3d histogram will be projected on the pads given by the following bounds
- // for rows and pads
- Int_t rowLow = iRow - rowRadius;
- UInt_t rowUp = iRow + rowRadius;
- Int_t padLow = iPad - padRadius;
- UInt_t padUp = iPad + padRadius;
- // if window goes out of chamber
- if (rowLow < 0) rowLow = 0;
- if (rowUp >= nRows) rowUp = nRows - 1;
- if (padLow < 0) padLow = 0;
- if (padUp >= nPads) padUp = nPads - 1;
-
- // project the histogram
- //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
- TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
-
- // get the number of entries in the spectrum
- Double_t entries = projH->GetEntries();
- if (entries < minEntries) { delete projH; continue; }
-
- // get the two calibration parameters mean of spectrum and RMS of spectrum
- Double_t histMean = projH->GetMean();
- Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
-
- // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
- Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
- Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
- if (fitResult != 0) {
- Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i.", chamber, rowLow, rowUp, padLow, padUp);
- delete projH;
- continue;
- }
-
- // get the two calibration parameters mean of gauss fit and sigma of gauss fit
- TF1* gausFit = projH->GetFunction("gaus");
- Double_t fitMean = gausFit->GetParameter(1);
- Double_t fitRMS = gausFit->GetParameter(2);
- Int_t numberFitPoints = gausFit->GetNumberFitPoints();
- if (numberFitPoints == 0) continue;
- Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
- delete gausFit;
- delete projH;
- if (fitMean <= 0) continue;
- //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
-
- // write the calibration parameters for each pad that the 3d histogram was projected onto
- // (with considering the step size) to the CalPads
- // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
- // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
- for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
- if (r < 0 || r >= (Int_t)nRows) continue;
- UInt_t nPads = roc.GetNPads(r);
- for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
- if (p < 0 || p >= (Int_t)nPads) continue;
- spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
- spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
- fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
- fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
- fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
- entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
-
- (*debugStream) << "calibKr" <<
- "sector=" << chamber << // chamber number
- "row=" << r << // row number
- "pad=" << p << // pad number
- "histMean=" << histMean << // mean of the spectrum
- "histRMS=" << histRMS << // RMS of the spectrum divided by the mean
- "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak
- "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak
- "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit
- "entries=" << entries << // number of entries for the spectrum
- "\n";
- }
- }
- }
- }
- }
-
- TFile f("calibKr.root", "recreate");
- spectrMeanCalPad->Write();
- spectrRMSCalPad->Write();
- fitMeanCalPad->Write();
- fitRMSCalPad->Write();
- fitNormChi2CalPad->Write();
- entriesCalPad->Write();
- f.Close();
- delete spectrMeanCalPad;
- delete spectrRMSCalPad;
- delete fitMeanCalPad;
- delete fitRMSCalPad;
- delete fitNormChi2CalPad;
- delete entriesCalPad;
- delete debugStream;
-}
-
-//_____________________________________________________________________
-TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
-{
- // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
- // replaces TH3F::ProjectZ() to gain more speed
-
- TAxis* xAxis = histo3D->GetXaxis();
- TAxis* yAxis = histo3D->GetYaxis();
- TAxis* zAxis = histo3D->GetZaxis();
- Double_t zMinVal = zAxis->GetXmin();
- Double_t zMaxVal = zAxis->GetXmax();
-
- Int_t nBinsZ = zAxis->GetNbins();
- TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
-
- Int_t nx = xAxis->GetNbins()+2;
- Int_t ny = yAxis->GetNbins()+2;
- Int_t bin = 0;
- Double_t entries = 0.;
- for (Int_t x = xMin; x <= xMax; x++) {
- for (Int_t y = yMin; y <= yMax; y++) {
- for (Int_t z = 0; z <= nBinsZ+1; z++) {
- bin = x + nx * (y + ny * z);
- Double_t val = histo3D->GetBinContent(bin);
- projH->Fill(zAxis->GetBinCenter(z), val);
- entries += val;
- }
- }
- }
- projH->SetEntries((Long64_t)entries);
- return projH;
-}
-
-//_____________________________________________________________________
-Long64_t AliTPCCalibKr::Merge(TCollection* list) {
-// merge function
-//
-cout << "Merge " << endl;
-
-if (!list)
-return 0;
-
-if (list->IsEmpty())
-return 1;
-
-TIterator* iter = list->MakeIterator();
-TObject* obj = 0;
-
- iter->Reset();
- Int_t count=0;
- while((obj = iter->Next()) != 0)
- {
- AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
- if (entry == 0) continue;
-
- for(int i=0; i<72; ++i) {
- if(IsCSide(i) == kTRUE && bCSide == kTRUE) {
- ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
- }
-
- if(IsCSide(i) == kFALSE && bASide == kTRUE) {
- ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
- }
- }
-
- count++;
- }
-
-return count;
-}
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes is hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+\r
+\r
+//Root includes\r
+#include <TH1F.h>\r
+#include <TH1D.h>\r
+#include <TH2F.h>\r
+#include <TH3F.h>\r
+#include <TString.h>\r
+#include <TMath.h>\r
+#include <TF1.h>\r
+#include <TRandom.h>\r
+#include <TDirectory.h>\r
+#include <TFile.h>\r
+#include <TAxis.h>\r
+//AliRoot includes\r
+#include "AliRawReader.h"\r
+#include "AliRawReaderRoot.h"\r
+#include "AliRawReaderDate.h"\r
+#include "AliTPCRawStream.h"\r
+#include "AliTPCCalROC.h"\r
+#include "AliTPCCalPad.h"\r
+#include "AliTPCROC.h"\r
+#include "AliMathBase.h"\r
+#include "TTreeStream.h"\r
+#include "AliTPCRawStreamFast.h"\r
+\r
+//date\r
+#include "event.h"\r
+\r
+//header file\r
+#include "AliTPCCalibKr.h"\r
+\r
+//----------------------------------------------------------------------------\r
+// The AliTPCCalibKr class description (TPC Kr calibration).\r
+//\r
+//\r
+// The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),\r
+// its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()). \r
+// \r
+// The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.\r
+// These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration \r
+// parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.\r
+// In addition the debugCalibKr.root file with debug information is created.\r
+//\r
+\r
+/*\r
+ \r
+// Usage example:\r
+//\r
+\r
+// 1. Analyse output histograms:\r
+TFile f("outHistFile.root");\r
+AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr")\r
+obj->Analyse();\r
+\r
+// 2. See calibration parameters e.g.:\r
+TFile f("calibKr.root");\r
+spectrMean->GetCalROC(70)->GetValue(40,40);\r
+fitMean->GetCalROC(70)->GetValue(40,40);\r
+\r
+// 3. See debug information e.g.:\r
+TFile f("debugCalibKr.root");\r
+.ls;\r
+\r
+// -- Print calibKr TTree content \r
+calibKr->Print();\r
+\r
+// -- Draw calibKr TTree variables\r
+calibKr.Draw("fitMean");\r
+\r
+*/\r
+\r
+\r
+//\r
+// Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)\r
+//-----------------------------------------------------------------------------\r
+\r
+ClassImp(AliTPCCalibKr)\r
+\r
+AliTPCCalibKr::AliTPCCalibKr() : \r
+ TObject(),\r
+ fASide(kTRUE),\r
+ fCSide(kTRUE),\r
+ fHistoKrArray(72),\r
+ fADCOverClustSizeMin(0.0),\r
+ fADCOverClustSizeMax(1.0e9),\r
+ fMaxADCOverClustADCMin(0.0),\r
+ fMaxADCOverClustADCMax(1.0e9),\r
+ fTimeMin(0.0),\r
+ fTimeMax(1.0e9),\r
+ fClustSizeMin(0.0),\r
+ fClustSizeMax(1.0e9),\r
+ fTimebinRmsIrocMin(0.0),\r
+ fPadRmsIrocMin(0.0),\r
+ fRowRmsIrocMin(0.0),\r
+ fClusterPadSize1DIrocMax(200),\r
+ fCurveCoefficientIroc(1.0e9),\r
+ fTimebinRmsOrocMin(0.0),\r
+ fPadRmsOrocMin(0.0),\r
+ fRowRmsOrocMin(0.0),\r
+ fClusterPadSize1DOrocMax(200),\r
+ fCurveCoefficientOroc(1.0e9)\r
+{\r
+ //\r
+ // default constructor\r
+ //\r
+\r
+ // TObjArray with histograms\r
+ fHistoKrArray.SetOwner(kTRUE); // is owner of histograms\r
+ fHistoKrArray.Clear();\r
+ \r
+ // init cuts (by Stefan)\r
+// SetADCOverClustSizeRange(7,200);\r
+// SetMaxADCOverClustADCRange(0.01,0.4);\r
+// SetTimeRange(200,600);\r
+// SetClustSizeRange(6,200);\r
+\r
+ //init cuts (by Adam)\r
+ SetTimebinRmsMin(1.6,0.8);\r
+ SetPadRmsMin(0.825,0.55);\r
+ SetRowRmsMin(0.1,0.1);\r
+ SetClusterPadSize1DMax(15,11);\r
+ SetCurveCoefficient(1.,2.);\r
+ \r
+ // init histograms\r
+ Init();\r
+}\r
+\r
+//_____________________________________________________________________\r
+AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : \r
+ TObject(pad),\r
+ \r
+ fASide(pad.fASide),\r
+ fCSide(pad.fCSide),\r
+ fHistoKrArray(72),\r
+ fADCOverClustSizeMin(pad.fADCOverClustSizeMin),\r
+ fADCOverClustSizeMax(pad.fADCOverClustSizeMax),\r
+ fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),\r
+ fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),\r
+ fTimeMin(pad.fTimeMin),\r
+ fTimeMax(pad.fTimeMax),\r
+ fClustSizeMin(pad.fClustSizeMin),\r
+ fClustSizeMax(pad.fClustSizeMax),\r
+ fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),\r
+ fPadRmsIrocMin(pad.fPadRmsIrocMin),\r
+ fRowRmsIrocMin(pad.fRowRmsIrocMin),\r
+ fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),\r
+ fCurveCoefficientIroc(pad.fCurveCoefficientIroc),\r
+ fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),\r
+ fPadRmsOrocMin(pad.fPadRmsOrocMin),\r
+ fRowRmsOrocMin(pad.fRowRmsOrocMin),\r
+ fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),\r
+ fCurveCoefficientOroc(pad.fCurveCoefficientOroc)\r
+{\r
+ // copy constructor\r
+ \r
+ // TObjArray with histograms\r
+ fHistoKrArray.SetOwner(kTRUE); // is owner of histograms\r
+ fHistoKrArray.Clear();\r
+\r
+ for (Int_t iSec = 0; iSec < 72; ++iSec) \r
+ {\r
+ TH3F *hOld = pad.GetHistoKr(iSec);\r
+ if(hOld) {\r
+ TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); \r
+ fHistoKrArray.AddAt(hNew,iSec);\r
+ }\r
+ }\r
+}\r
+\r
+//_____________________________________________________________________\r
+AliTPCCalibKr::~AliTPCCalibKr() \r
+{\r
+ //\r
+ // destructor\r
+ //\r
+\r
+ // for (Int_t iSec = 0; iSec < 72; ++iSec) {\r
+ // if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);\r
+ // }\r
+ fHistoKrArray.Delete();\r
+}\r
+\r
+//_____________________________________________________________________\r
+AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)\r
+{\r
+ // assignment operator\r
+\r
+ if (&source == this) return *this;\r
+ new (this) AliTPCCalibKr(source);\r
+\r
+ return *this;\r
+}\r
+\r
+//_____________________________________________________________________\r
+void AliTPCCalibKr::Init()\r
+{\r
+ // \r
+ // init output histograms \r
+ //\r
+\r
+ // add histograms to the TObjArray\r
+ for(Int_t i=0; i<72; ++i) {\r
+ \r
+ // C - side\r
+ if(IsCSide(i) == kTRUE && fCSide == kTRUE) {\r
+ TH3F *hist = CreateHisto(i);\r
+ if(hist) fHistoKrArray.AddAt(hist,i);\r
+ }\r
+ \r
+ // A - side\r
+ if(IsCSide(i) == kFALSE && fASide == kTRUE) {\r
+ TH3F *hist = CreateHisto(i);\r
+ if(hist) fHistoKrArray.AddAt(hist,i);\r
+ }\r
+ }\r
+}\r
+ \r
+//_____________________________________________________________________\r
+Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)\r
+{\r
+ //\r
+ // process events \r
+ // call event by event\r
+ //\r
+\r
+ if(cluster) Update(cluster);\r
+ else return kFALSE;\r
+ \r
+ return kTRUE;\r
+}\r
+\r
+//_____________________________________________________________________\r
+TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)\r
+{\r
+ //\r
+ // create new histogram\r
+ //\r
+ char name[256];\r
+ TH3F *h;\r
+\r
+ sprintf(name,"ADCcluster_ch%d",chamber);\r
+\r
+ if( IsIROC(chamber) == kTRUE ) \r
+ {\r
+ h = new TH3F(name,name,63,0,63,108,0,108,200,100,2500);\r
+ } else {\r
+ h = new TH3F(name,name,96,0,96,140,0,140,200,100,1700);\r
+ }\r
+ h->SetXTitle("padrow");\r
+ h->SetYTitle("pad");\r
+ h->SetZTitle("fADC");\r
+\r
+return h;\r
+}\r
+\r
+//_____________________________________________________________________\r
+Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)\r
+{\r
+// check if IROCs\r
+// returns kTRUE if IROCs and kFALSE if OROCs \r
+\r
+ if(chamber>=0 && chamber<36) return kTRUE;\r
+\r
+return kFALSE;\r
+}\r
+\r
+//_____________________________________________________________________\r
+Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)\r
+{\r
+// check if C side\r
+// returns kTRUE if C side and kFALSE if A side\r
+\r
+ if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;\r
+\r
+return kFALSE;\r
+}\r
+ \r
+//_____________________________________________________________________\r
+Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)\r
+{\r
+ //\r
+ // fill existing histograms\r
+ //\r
+\r
+ if (!Accept(cl)) return kFALSE;\r
+ TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());\r
+ if(!h) return kFALSE;\r
+ \r
+ h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());\r
+ \r
+ return kTRUE;\r
+}\r
+\r
+//_____________________________________________________________________\r
+Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){\r
+ //\r
+ // cuts\r
+ //\r
+ /*\r
+ TCut cutR0("cutR0","fADCcluster/fSize<200"); // adjust it according v seetings - \r
+ TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal\r
+ TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4"); // digital noise removal\r
+ TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal\r
+ TCut cutR4("cutR4","fMax.fTime>200"); // noise removal\r
+ TCut cutR5("cutR5","fMax.fTime<600"); // noise removal\r
+ TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks\r
+ TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;\r
+ */\r
+ /*\r
+ //R0\r
+ if ((float)cl->GetADCcluster()/ cl->GetSize() >200) return kFALSE;\r
+ // R1\r
+ if ((float)cl->GetADCcluster()/ cl->GetSize() <7) return kFALSE;\r
+ //R2\r
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4) return kFALSE;\r
+ //R3\r
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;\r
+ //R4\r
+ if (cl->GetMax().GetTime() < 200) return kFALSE;\r
+ //R5\r
+ if (cl->GetMax().GetTime() > 600) return kFALSE;\r
+ //S1\r
+ if (cl->GetSize()>200) return kFALSE;\r
+ if (cl->GetSize()<6) return kFALSE;\r
+\r
+ SetADCOverClustSizeRange(7,200);\r
+ SetMaxADCOverClustADCRange(0.01,0.4);\r
+ SetTimeRange(200,600);\r
+ SetClustSizeRange(6,200);\r
+ */\r
+\r
+ //R0\r
+ if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax) return kFALSE;\r
+ // R1\r
+ if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin) return kFALSE;\r
+ //R2\r
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax) return kFALSE;\r
+ //R3\r
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;\r
+ //R4\r
+ if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;\r
+ //R5\r
+ if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;\r
+ //S1\r
+ if (cl->GetSize()>fClustSizeMax) return kFALSE;\r
+ if (cl->GetSize()<fClustSizeMin) return kFALSE;\r
+\r
+ //\r
+ // cuts by Adam\r
+ //\r
+ /*\r
+ TCut cutAI0("cutAI0","fTimebinRMS>1.6");\r
+ TCut cutAI1("cutAI1","fPadRMS>0.825"); \r
+ TCut cutAI2("cutAI2","fRowRMS>0.1"); \r
+ TCut cutAI3("cutAI3","fPads1D<15"); \r
+ TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0"); \r
+\r
+ TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;\r
+\r
+ TCut cutAO0("cutAO0","fTimebinRMS>0.8");\r
+ TCut cutAO1("cutAO1","fPadRMS>0.55"); \r
+ TCut cutAO2("cutAO2","fRowRMS>0.1"); \r
+ TCut cutAO3("cutAO3","fPads1D<11"); \r
+ TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0"); \r
+\r
+ TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;\r
+ TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");\r
+\r
+ */\r
+ if(cl->GetSec()<36){ //IROCs\r
+ //AI0\r
+ if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;\r
+ //AI1\r
+ if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;\r
+ //AI2\r
+ if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;\r
+ //AI3\r
+ if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;\r
+ //AI4\r
+ if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;\r
+ }else{ //OROCs\r
+ //AO0\r
+ if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;\r
+ //AO1\r
+ if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;\r
+ //AO2\r
+ if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;\r
+ //AO3\r
+ if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;\r
+ //AO4\r
+ if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;\r
+ }\r
+\r
+ return kTRUE;\r
+\r
+}\r
+\r
+//_____________________________________________________________________\r
+TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const\r
+{\r
+ // get histograms from fHistoKrArray\r
+ return (TH3F*) fHistoKrArray.At(chamber);\r
+}\r
+ \r
+//_____________________________________________________________________\r
+void AliTPCCalibKr::Analyse() \r
+{\r
+ //\r
+ // analyse the histograms and extract krypton calibration parameters\r
+ //\r
+\r
+ // AliTPCCalPads that will contain the calibration parameters\r
+ AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");\r
+ AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");\r
+ AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");\r
+ AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");\r
+ AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");\r
+ AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");\r
+\r
+ // file stream for debugging purposes\r
+ TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");\r
+\r
+ // if entries in spectrum less than minEntries, then the fit won't be performed\r
+ Int_t minEntries = 1; //300;\r
+\r
+ Double_t windowFrac = 0.12;\r
+ // the 3d histogram will be projected on the pads given by the following window size\r
+ // set the numbers to 0 if you want to do a pad-by-pad calibration\r
+ UInt_t rowRadius = 2;\r
+ UInt_t padRadius = 4;\r
+ \r
+ // the step size by which pad and row are incremented is given by the following numbers\r
+ // set them to 1 if you want the finest granularity\r
+ UInt_t rowStep = 1; // formerly: 2*rowRadius\r
+ UInt_t padStep = 1; // formerly: 2*padRadius\r
+\r
+ for (Int_t chamber = 0; chamber < 72; chamber++) {\r
+ //if (chamber != 71) continue;\r
+ AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()\r
+ \r
+ // Usually I would traverse each pad, take the spectrum for its neighbourhood and\r
+ // obtain the calibration parameters. This takes very long, so instead I assign the same\r
+ // calibration values to the whole neighbourhood and then go on to the next neighbourhood.\r
+ UInt_t nRows = roc.GetNrows();\r
+ for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {\r
+ UInt_t nPads = roc.GetNPads(iRow);\r
+ //if (iRow >= 10) break;\r
+ for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {\r
+ //if (iPad >= 20) break;\r
+ TH3F* h = GetHistoKr(chamber);\r
+ if (!h) continue;\r
+ \r
+ // the 3d histogram will be projected on the pads given by the following bounds\r
+ // for rows and pads\r
+ Int_t rowLow = iRow - rowRadius;\r
+ UInt_t rowUp = iRow + rowRadius;\r
+ Int_t padLow = iPad - padRadius;\r
+ UInt_t padUp = iPad + padRadius;\r
+ // if window goes out of chamber\r
+ if (rowLow < 0) rowLow = 0;\r
+ if (rowUp >= nRows) rowUp = nRows - 1;\r
+ if (padLow < 0) padLow = 0;\r
+ if (padUp >= nPads) padUp = nPads - 1;\r
+\r
+ // project the histogram\r
+ //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW\r
+ TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);\r
+ \r
+ // get the number of entries in the spectrum\r
+ Double_t entries = projH->GetEntries();\r
+ if (entries < minEntries) { delete projH; continue; }\r
+ \r
+ // get the two calibration parameters mean of spectrum and RMS of spectrum\r
+ Double_t histMean = projH->GetMean();\r
+ Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;\r
+ \r
+ // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted\r
+ Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());\r
+ Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);\r
+ Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);\r
+ Double_t integCharge = projH->Integral(minBin,maxBin); \r
+\r
+ Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);\r
+\r
+ if (fitResult != 0) {\r
+ Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);\r
+ //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, %f\n", chamber, iRow, iPad, entries, maxEntries);\r
+ \r
+ delete projH;\r
+ continue;\r
+ }\r
+ \r
+ // get the two calibration parameters mean of gauss fit and sigma of gauss fit\r
+ TF1* gausFit = projH->GetFunction("gaus");\r
+ Double_t fitMean = gausFit->GetParameter(1);\r
+ Double_t fitRMS = gausFit->GetParameter(2);\r
+ Int_t numberFitPoints = gausFit->GetNumberFitPoints();\r
+ if (numberFitPoints == 0) continue;\r
+ Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;\r
+ delete gausFit;\r
+ delete projH;\r
+ if (fitMean <= 0) continue;\r
+ //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);\r
+ \r
+ // write the calibration parameters for each pad that the 3d histogram was projected onto\r
+ // (with considering the step size) to the CalPads\r
+ // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions\r
+ // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction\r
+ for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {\r
+ if (r < 0 || r >= (Int_t)nRows) continue;\r
+ UInt_t nPadsR = roc.GetNPads(r);\r
+ for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {\r
+ if (p < 0 || p >= (Int_t)nPadsR) continue;\r
+ spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);\r
+ spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);\r
+ fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);\r
+ fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);\r
+ fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);\r
+ entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);\r
+\r
+ (*debugStream) << "calibKr" <<\r
+ "sector=" << chamber << // chamber number\r
+ "row=" << r << // row number\r
+ "pad=" << p << // pad number\r
+ "histMean=" << histMean << // mean of the spectrum\r
+ "histRMS=" << histRMS << // RMS of the spectrum divided by the mean\r
+ "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak\r
+ "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak\r
+ "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit\r
+ "entries=" << entries << // number of entries for the spectrum\r
+ "\n";\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ TFile f("calibKr.root", "recreate");\r
+ spectrMeanCalPad->Write();\r
+ spectrRMSCalPad->Write();\r
+ fitMeanCalPad->Write();\r
+ fitRMSCalPad->Write();\r
+ fitNormChi2CalPad->Write();\r
+ entriesCalPad->Write();\r
+ f.Close();\r
+ delete spectrMeanCalPad;\r
+ delete spectrRMSCalPad;\r
+ delete fitMeanCalPad;\r
+ delete fitRMSCalPad;\r
+ delete fitNormChi2CalPad;\r
+ delete entriesCalPad;\r
+ delete debugStream;\r
+}\r
+\r
+//_____________________________________________________________________\r
+TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)\r
+{\r
+ // project the z-axis of a 3d histo to a specific range of the x- and y-axes,\r
+ // replaces TH3F::ProjectZ() to gain more speed\r
+\r
+ TAxis* xAxis = histo3D->GetXaxis();\r
+ TAxis* yAxis = histo3D->GetYaxis();\r
+ TAxis* zAxis = histo3D->GetZaxis();\r
+ Double_t zMinVal = zAxis->GetXmin();\r
+ Double_t zMaxVal = zAxis->GetXmax();\r
+ \r
+ Int_t nBinsZ = zAxis->GetNbins();\r
+ TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);\r
+\r
+ Int_t nx = xAxis->GetNbins()+2;\r
+ Int_t ny = yAxis->GetNbins()+2;\r
+ Int_t bin = 0;\r
+ Double_t entries = 0.;\r
+ for (Int_t x = xMin; x <= xMax; x++) {\r
+ for (Int_t y = yMin; y <= yMax; y++) {\r
+ for (Int_t z = 0; z <= nBinsZ+1; z++) {\r
+ bin = x + nx * (y + ny * z);\r
+ Double_t val = histo3D->GetBinContent(bin);\r
+ projH->Fill(zAxis->GetBinCenter(z), val);\r
+ entries += val;\r
+ }\r
+ }\r
+ }\r
+ projH->SetEntries((Long64_t)entries);\r
+ return projH;\r
+}\r
+\r
+//_____________________________________________________________________\r
+Long64_t AliTPCCalibKr::Merge(TCollection* list) {\r
+// merge function \r
+//\r
+cout << "Merge " << endl;\r
+\r
+if (!list)\r
+return 0;\r
+\r
+if (list->IsEmpty())\r
+return 1;\r
+\r
+TIterator* iter = list->MakeIterator();\r
+TObject* obj = 0;\r
+\r
+ iter->Reset();\r
+ Int_t count=0;\r
+ while((obj = iter->Next()) != 0)\r
+ {\r
+ AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);\r
+ if (entry == 0) continue; \r
+\r
+ for(int i=0; i<72; ++i) { \r
+ if(IsCSide(i) == kTRUE && fCSide == kTRUE) { \r
+ ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) ); \r
+ } \r
+\r
+ if(IsCSide(i) == kFALSE && fASide == kTRUE) { \r
+ ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) ); \r
+ } \r
+ } \r
+\r
+ count++;\r
+ }\r
+\r
+return count;\r
+}\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff