Gain calibration per chamber:
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Apr 2012 13:14:13 +0000 (13:14 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Apr 2012 13:14:13 +0000 (13:14 +0000)
+++ AliTPCcalibGainMult.h                 - +  TGraphErrors* GetGainPerChamber(Int_t padRegion=1, Bool_t plotQA=kFALSE);
+++ AliTPCcalibGainMult.cxx               - +  Export of the gain per chamber + Enabling  ProcessV0 and ProcessTOF for identified particles selection
+++ AliTPCseed.cxx      (working copy)    - +  Use of gain
+++ AliTPCPreprocessorOffline.h (working copy) + AnalyzeGainChamberByChamber();
+++ AliTPCPreprocessorOffline.cxx       (working copy)  AnalyzeGainChamberByChamber();
+

TPC/AliTPCPreprocessorOffline.cxx
TPC/AliTPCPreprocessorOffline.h
TPC/AliTPCcalibGainMult.cxx
TPC/AliTPCcalibGainMult.h
TPC/AliTPCseed.cxx

index c79ab80..41185d4 100644 (file)
@@ -838,6 +838,7 @@ void AliTPCPreprocessorOffline::CalibTimeGain(const Char_t* fileName, Int_t star
   AnalyzeAttachment(startRunNumber,endRunNumber);
   AnalyzePadRegionGain();
   AnalyzeGainMultiplicity();
+  AnalyzeGainChamberByChamber();
   //
   // 3. Make control plots
   //
@@ -1143,6 +1144,26 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
 
 }
 
+Bool_t AliTPCPreprocessorOffline::AnalyzeGainChamberByChamber(){
+  //
+  // get chamber by chamber gain
+  //
+  TGraphErrors *grShort  = fGainMult->GetGainPerChamber(0);
+  TGraphErrors *grMedium = fGainMult->GetGainPerChamber(1);
+  TGraphErrors *grLong   = fGainMult->GetGainPerChamber(2);
+  if (grShort==0x0 || grMedium==0x0 || grLong==0x0) {
+    delete grShort;
+    delete grMedium;
+    delete grLong;
+    return kFALSE;
+  }
+
+  fGainArray->AddLast(grShort);
+  fGainArray->AddLast(grMedium);
+  fGainArray->AddLast(grLong);
+
+  return kTRUE;
+}
 
 void AliTPCPreprocessorOffline::UpdateOCDBGain(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath){
   //
index 22e3dd2..5613080 100644 (file)
@@ -50,6 +50,7 @@ public:
   Bool_t AnalyzeAttachment(Int_t startRunNumber, Int_t endRunNumber, Int_t minEntriesFit = 2000);
   Bool_t AnalyzePadRegionGain();
   Bool_t AnalyzeGainMultiplicity();
+  Bool_t AnalyzeGainChamberByChamber();
   void SetTimeGainRange(Double_t minGain=2.0, Double_t maxGain = 3.0) 
        {fMinGain = minGain; fMaxGain = maxGain;};
   Bool_t ValidateTimeGain();
index cd048a3..f12e22b 100644 (file)
@@ -25,6 +25,7 @@ Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
 
 
 #include "Riostream.h"
+#include "TROOT.h"
 #include "TChain.h"
 #include "TTree.h"
 #include "TH1F.h"
@@ -36,6 +37,7 @@ Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
 #include "TF1.h"
 #include "TVectorF.h"
 #include "TProfile.h"
+#include "TGraphErrors.h"
 
 #include "AliTPCcalibDB.h"
 #include "AliTPCclusterMI.h"
@@ -65,7 +67,6 @@ Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
 #include "AliTracker.h"
 #include "AliTPCTransform.h"
 #include "AliTPCROC.h"
-#include "TROOT.h"
 
 ClassImp(AliTPCcalibGainMult)
 
@@ -276,8 +277,8 @@ void AliTPCcalibGainMult::Process(AliESDEvent *event) {
     return;
   }
   //
-  //ProcessV0s(event);   // 
-  //ProcessTOF(event);   //
+  ProcessV0s(event);   // 
+  ProcessTOF(event);   //
   //ProcessKinks(event); // not relyable
   DumpHPT(event);      // 
   UInt_t runNumber = event->GetRunNumber();
@@ -991,7 +992,7 @@ void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftra
   if (isKaon||isProton||isPion||isElectron||isMIP||isMuon) isSelected=kTRUE;
   isHighPt = kFALSE;
   if (!isSelected) isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005);  
-  if (counter>kMax && ((1/track->Pt()*gRandom->Rndm())>kMax/counter)) return; 
+  //if (counter>kMax && ((1/track->Pt()*gRandom->Rndm())>kMax/counter)) return; 
   isSelected|=isHighPt;
   //
   //
@@ -1083,10 +1084,10 @@ void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftra
       "pout="<<pout<<
       "tgl="<<tgl<<
       "track.="<<track<<
-      "trackIn.="<<trackIn<<
-      "trackOut.="<<trackOut<<
-      "tpcOut.="<<tpcOut<<
-      "seed.="<<seed<<
+      //"trackIn.="<<trackIn<<
+      //"trackOut.="<<trackOut<<
+      //"tpcOut.="<<tpcOut<<
+      //"seed.="<<seed<<
       "medianMIP0="<<medianMIP0<<    // median MIP position as estimated from the array of (not only) "MIPS"
       //dedx 
       "truncUp.="<<&truncUp<<
@@ -1781,6 +1782,103 @@ void AliTPCcalibGainMult::ProcessTOF(const AliESDEvent * event){
 }
 
 
+TGraphErrors* AliTPCcalibGainMult::GetGainPerChamber(Int_t padRegion/*=1*/, Bool_t plotQA/*=kFALSE*/)
+{
+  //
+  // Extract gain variations per chamger for 'padRegion'
+  //
+
+  if (padRegion<0||padRegion>2) return 0x0;
+  
+  fHistGainSector->GetAxis(2)->SetRangeUser(padRegion,padRegion);
+  TH2D * histGainSec = fHistGainSector->Projection(0,1);
+//   TH1D* proj=fHistGainSector->Projection(0);
+//   Double_t max=proj->GetBinCenter(proj->GetMaximumBin());
+//   delete proj;
+  //
+  TObjArray arr;
+//   TF1 fg("gaus","gaus",histGainSec->GetYaxis()->GetXmin()+1,histGainSec->GetYaxis()->GetXmax()-1);
+//   histGainSec->FitSlicesY(&fg, 0, -1, 0, "QNR", &arr);
+  histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
+  TH1D * meanGainSec = (TH1D*)arr.At(1);
+  Double_t gainsIROC[36]={0.};
+  Double_t gainsOROC[36]={0.};
+  Double_t gains[72]={0.};
+  Double_t gainsErr[72]={0.};
+  TGraphErrors *gr=new TGraphErrors(36);
+  //
+  for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
+    TH1D *slice=histGainSec->ProjectionY("_py",isec,isec);
+    Double_t max=slice->GetBinCenter(slice->GetMaximumBin());
+    TF1 fg("gaus","gaus",max-30,max+30);
+    slice->Fit(&fg,"QNR");
+    meanGainSec->SetBinContent(isec,fg.GetParameter(1));
+    meanGainSec->SetBinError(isec,fg.GetParError(1));
+    
+//     cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
+    gains[isec-1] = meanGainSec->GetBinContent(isec);
+    if (isec < 37) {
+      gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
+    } else {
+      gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
+    }
+    gainsErr[isec-1]=meanGainSec->GetBinError(isec);
+    delete slice;
+  }
+  Double_t meanIroc = TMath::Median(36, gainsIROC);
+  Double_t meanOroc = TMath::Median(36, gainsOROC);
+  if (TMath::Abs(meanIroc)<1e-30) meanIroc=1.;
+  if (TMath::Abs(meanOroc)<1e-30) meanOroc=1.;
+  for(Int_t i = 0; i < 36; i++) {
+    gains[i] /= meanIroc;
+    gainsErr[i] /= meanIroc;
+  }
+  
+  for(Int_t i = 36; i < 72; i++){
+    gains[i] /= meanOroc;
+    gainsErr[i] /= meanOroc;
+  }
+  //
+  Int_t ipoint=0;
+  for(Int_t i = 0; i < 72; i++) {
+    if (padRegion==0 && i>35) continue;
+    if ( (padRegion==1 || padRegion==2) && i<36) continue;
+
+    if (gains[i]<1e-20 || gainsErr[i]/gains[i]>.2){
+      AliWarning(Form("Invalid chamber gain in ROC/region: %d / %d", i, padRegion));
+      gains[i]=1;
+      gainsErr[i]=1;
+    }
+    
+    gr->SetPoint(ipoint,i,gains[i]);
+    gr->SetPointError(ipoint,0,gainsErr[i]);
+    ++ipoint;
+  }
+
+  const char* names[3]={"SHORT","MEDIUM","LONG"};
+  gr->SetNameTitle(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]),Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]));
+
+
+  //=====================================
+  // Do QA plotting if requested
+  if (plotQA){
+    TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("cQA");
+    if (!c) c=new TCanvas("cQA","cQA");
+    c->Clear();
+    c->Divide(1,2);
+    c->cd(1);
+    histGainSec->DrawCopy("colz");
+    meanGainSec->DrawCopy("same");
+    gr->SetMarkerStyle(20);
+    gr->SetMarkerSize(.5);
+    c->cd(2);
+    gr->Draw("alp");
+  }
+  
+  delete histGainSec;
+  return gr;
+}
+
 // void AliTPCcalibGainMult::Terminate(){
 //   //
 //   // Terminate function
index cd56d81..2ad8d39 100644 (file)
@@ -50,6 +50,7 @@ public:
   THnSparseF * GetHistdEdxTot() const { return fHistdEdxTot;}      // 4D dedx histogram
   TTree *      GetdEdxTree() const {return fdEdxTree;}         // tree for the later minimization
 
+  TGraphErrors* GetGainPerChamber(Int_t padRegion=1, Bool_t plotQA=kFALSE);
   //
   void SetMIPvalue(Float_t mip){fMIP = mip;};
   void SetLowerTrunc(Float_t lowerTrunc){fLowerTrunc = lowerTrunc;};
index 4dd8143..9bbe38f 100644 (file)
@@ -1143,6 +1143,7 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   Float_t corrTimeGain = 1;
   TObjArray * timeGainSplines = 0x0;
   TGraphErrors * grPadEqual = 0x0;
+  TGraphErrors*  grChamberGain[3]={0x0,0x0,0x0};
   //
   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
@@ -1165,8 +1166,11 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
        //
        if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
        if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
+        const char* names[3]={"SHORT","MEDIUM","LONG"};
+        for (Int_t iPadRegion=0; iPadRegion<3; ++iPadRegion)
+          grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion]));
       }
-  }   
+  }
   
   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
@@ -1258,13 +1262,19 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
     // pad region equalization outside of cluster param
     //
     Float_t gainEqualPadRegion = 1;
-    if (grPadEqual) gainEqualPadRegion = grPadEqual->Eval(ipad);
+    if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad);
+    //
+    // chamber-by-chamber equalization outside gain map
+    //
+    Float_t gainChamber = 1;
+    if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector());
     //
     amp[ncl]=charge;
     amp[ncl]/=gainGG;
     amp[ncl]/=gainPad;
     amp[ncl]/=corrPos;
     amp[ncl]/=gainEqualPadRegion;
+    amp[ncl]/=gainChamber;
     //
     ncl++;
   }