-/**************************************************************************\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 "AliTPCCalROC.h"\r
-#include "AliTPCCalPad.h"\r
-#include "AliTPCROC.h"\r
-#include "AliMathBase.h"\r
-#include "TTreeStream.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->SetRadius(0,0);\r
-obj->SetStep(1,1);\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
- fIrocHistogramMin(100.),\r
- fIrocHistogramMax(6000.),\r
- fIrocHistogramNbins(200),\r
- fOrocHistogramMin(100.),\r
- fOrocHistogramMax(5500.),\r
- fOrocHistogramNbins(200),\r
- fRowRadius(0),\r
- fPadRadius(0),\r
- fRowStep(1),\r
- fPadStep(1)\r
-\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
- //set histograms settings\r
- SetIrocHistogram(200,100,6000);\r
- SetOrocHistogram(200,100,5500);\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
- fIrocHistogramMin(pad.fIrocHistogramMin),\r
- fIrocHistogramMax(pad.fIrocHistogramMax),\r
- fIrocHistogramNbins(pad.fIrocHistogramNbins),\r
- fOrocHistogramMin(pad.fOrocHistogramMin),\r
- fOrocHistogramMax(pad.fOrocHistogramMax),\r
- fOrocHistogramNbins(pad.fOrocHistogramNbins),\r
- fRowRadius(pad.fRowRadius),\r
- fPadRadius(pad.fPadRadius),\r
- fRowStep(pad.fRowStep),\r
- fPadStep(pad.fPadStep)\r
-\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
- snprintf(name,256,"ADCcluster_ch%d",chamber);\r
-\r
- if( IsIROC(chamber) == kTRUE ) \r
- {\r
- h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);\r
- } else {\r
- h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);\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 = fRowRadius;//4;\r
- UInt_t padRadius = fPadRadius;//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 = fRowStep;//1;//2 // formerly: 2*rowRadius\r
- UInt_t padStep = fPadStep;//1;//4 // 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 + rowStep-1;\r
- Int_t padLow = iPad - padRadius;\r
- UInt_t padUp = iPad + padRadius + padStep-1;\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+1, rowUp+1, padLow+1, padUp+1); // SLOW\r
- TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);\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
-\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
+/**************************************************************************
+ * 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 "AliTPCCalROC.h"
+#include "AliTPCCalPad.h"
+#include "AliTPCROC.h"
+#include "AliMathBase.h"
+#include "TTreeStream.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->SetRadius(0,0);
+obj->SetStep(1,1);
+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(),
+ fASide(kTRUE),
+ fCSide(kTRUE),
+ fHistoKrArray(72),
+ fADCOverClustSizeMin(0.0),
+ fADCOverClustSizeMax(1.0e9),
+ fMaxADCOverClustADCMin(0.0),
+ fMaxADCOverClustADCMax(1.0e9),
+ fTimeMin(0.0),
+ fTimeMax(1.0e9),
+ fClustSizeMin(0.0),
+ fClustSizeMax(1.0e9),
+ fTimebinRmsIrocMin(0.0),
+ fPadRmsIrocMin(0.0),
+ fRowRmsIrocMin(0.0),
+ fClusterPadSize1DIrocMax(200),
+ fCurveCoefficientIroc(1.0e9),
+ fTimebinRmsOrocMin(0.0),
+ fPadRmsOrocMin(0.0),
+ fRowRmsOrocMin(0.0),
+ fClusterPadSize1DOrocMax(200),
+ fCurveCoefficientOroc(1.0e9),
+ fIrocHistogramMin(100.),
+ fIrocHistogramMax(6000.),
+ fIrocHistogramNbins(200),
+ fOrocHistogramMin(100.),
+ fOrocHistogramMax(5500.),
+ fOrocHistogramNbins(200),
+ fRowRadius(0),
+ fPadRadius(0),
+ fRowStep(1),
+ fPadStep(1)
+
+{
+ //
+ // default constructor
+ //
+
+ // TObjArray with histograms
+ fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
+ fHistoKrArray.Clear();
+
+ // init cuts (by Stefan)
+// SetADCOverClustSizeRange(7,200);
+// SetMaxADCOverClustADCRange(0.01,0.4);
+// SetTimeRange(200,600);
+// SetClustSizeRange(6,200);
+
+ //init cuts (by Adam)
+ SetTimebinRmsMin(1.6,0.8);
+ SetPadRmsMin(0.825,0.55);
+ SetRowRmsMin(0.1,0.1);
+ SetClusterPadSize1DMax(15,11);
+ SetCurveCoefficient(1.,2.);
+
+ //set histograms settings
+ SetIrocHistogram(200,100,6000);
+ SetOrocHistogram(200,100,5500);
+
+ // init histograms
+ //Init();
+}
+
+//_____________________________________________________________________
+AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
+ TObject(pad),
+
+ fASide(pad.fASide),
+ fCSide(pad.fCSide),
+ fHistoKrArray(72),
+ fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
+ fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
+ fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
+ fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
+ fTimeMin(pad.fTimeMin),
+ fTimeMax(pad.fTimeMax),
+ fClustSizeMin(pad.fClustSizeMin),
+ fClustSizeMax(pad.fClustSizeMax),
+ fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),
+ fPadRmsIrocMin(pad.fPadRmsIrocMin),
+ fRowRmsIrocMin(pad.fRowRmsIrocMin),
+ fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),
+ fCurveCoefficientIroc(pad.fCurveCoefficientIroc),
+ fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),
+ fPadRmsOrocMin(pad.fPadRmsOrocMin),
+ fRowRmsOrocMin(pad.fRowRmsOrocMin),
+ fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),
+ fCurveCoefficientOroc(pad.fCurveCoefficientOroc),
+ fIrocHistogramMin(pad.fIrocHistogramMin),
+ fIrocHistogramMax(pad.fIrocHistogramMax),
+ fIrocHistogramNbins(pad.fIrocHistogramNbins),
+ fOrocHistogramMin(pad.fOrocHistogramMin),
+ fOrocHistogramMax(pad.fOrocHistogramMax),
+ fOrocHistogramNbins(pad.fOrocHistogramNbins),
+ fRowRadius(pad.fRowRadius),
+ fPadRadius(pad.fPadRadius),
+ fRowStep(pad.fRowStep),
+ fPadStep(pad.fPadStep)
+
+{
+ // copy constructor
+
+ // TObjArray with histograms
+ fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
+ fHistoKrArray.Clear();
+
+ 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) {
+ // if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
+ // }
+ 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 output histograms
+ //
+
+ // add histograms to the TObjArray
+ for(Int_t i=0; i<72; ++i) {
+
+ // C - side
+ if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
+ TH3F *hist = CreateHisto(i);
+ if(hist) fHistoKrArray.AddAt(hist,i);
+ }
+
+ // A - side
+ if(IsCSide(i) == kFALSE && fASide == 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;
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
+{
+ //
+ // create new histogram
+ //
+ char name[256];
+ TH3F *h;
+
+ snprintf(name,256,"ADCcluster_ch%d",chamber);
+
+ if( IsIROC(chamber) == kTRUE )
+ {
+ h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);
+ } else {
+ h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);
+ }
+ 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 cutR4("cutR4","fMax.fTime>200"); // noise removal
+ TCut cutR5("cutR5","fMax.fTime<600"); // noise removal
+ TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
+ TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
+ */
+ /*
+ //R0
+ if ((float)cl->GetADCcluster()/ cl->GetSize() >200) return kFALSE;
+ // R1
+ if ((float)cl->GetADCcluster()/ cl->GetSize() <7) return kFALSE;
+ //R2
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4) return kFALSE;
+ //R3
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
+ //R4
+ if (cl->GetMax().GetTime() < 200) return kFALSE;
+ //R5
+ if (cl->GetMax().GetTime() > 600) return kFALSE;
+ //S1
+ if (cl->GetSize()>200) return kFALSE;
+ if (cl->GetSize()<6) return kFALSE;
+
+ SetADCOverClustSizeRange(7,200);
+ SetMaxADCOverClustADCRange(0.01,0.4);
+ SetTimeRange(200,600);
+ SetClustSizeRange(6,200);
+ */
+
+ //R0
+ if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax) return kFALSE;
+ // R1
+ if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin) return kFALSE;
+ //R2
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax) return kFALSE;
+ //R3
+ if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
+ //R4
+ if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
+ //R5
+ if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
+ //S1
+ if (cl->GetSize()>fClustSizeMax) return kFALSE;
+ if (cl->GetSize()<fClustSizeMin) return kFALSE;
+
+ //
+ // cuts by Adam
+ //
+ /*
+ TCut cutAI0("cutAI0","fTimebinRMS>1.6");
+ TCut cutAI1("cutAI1","fPadRMS>0.825");
+ TCut cutAI2("cutAI2","fRowRMS>0.1");
+ TCut cutAI3("cutAI3","fPads1D<15");
+ TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0");
+
+ TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;
+
+ TCut cutAO0("cutAO0","fTimebinRMS>0.8");
+ TCut cutAO1("cutAO1","fPadRMS>0.55");
+ TCut cutAO2("cutAO2","fRowRMS>0.1");
+ TCut cutAO3("cutAO3","fPads1D<11");
+ TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0");
+
+ TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;
+ TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");
+
+ */
+ if(cl->GetSec()<36){ //IROCs
+ //AI0
+ if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;
+ //AI1
+ if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;
+ //AI2
+ if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;
+ //AI3
+ if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;
+ //AI4
+ if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
+ }else{ //OROCs
+ //AO0
+ if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;
+ //AO1
+ if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;
+ //AO2
+ if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;
+ //AO3
+ if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;
+ //AO4
+ if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) 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 = fRowRadius;//4;
+ UInt_t padRadius = fPadRadius;//4;
+
+ // 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 = fRowStep;//1;//2 // formerly: 2*rowRadius
+ UInt_t padStep = fPadStep;//1;//4 // 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 + rowStep-1;
+ Int_t padLow = iPad - padRadius;
+ UInt_t padUp = iPad + padRadius + padStep-1;
+ // 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+1, rowUp+1, padLow+1, padUp+1); // SLOW
+ TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);
+
+ // 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 minBin = projH->FindBin((1.-windowFrac) * maxEntries);
+ Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
+ Double_t integCharge = projH->Integral(minBin,maxBin);
+
+ 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, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
+ //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, %f\n", chamber, iRow, iPad, entries, maxEntries);
+
+ 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 nPadsR = 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)nPadsR) 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
+//
+
+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 && fCSide == kTRUE) {
+ ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
+ }
+
+ if(IsCSide(i) == kFALSE && fASide == kTRUE) {
+ ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
+ }
+ }
+
+ count++;
+ }
+
+return count;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff