]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibGainMult.cxx
fix compilation of results macro
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibGainMult.cxx
index 10e0223a21c6714d6459d15cfbc73729f8fb05aa..0ec6c5e12702b7f5d28882f797e49be51cfbe1e4 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)
 
@@ -147,14 +148,14 @@ AliTPCcalibGainMult::AliTPCcalibGainMult(const Text_t *name, const Text_t *title
   //
   //
   //
-  Int_t binsPadEqual[6]    = { 400, 400,    4,   20,   50, 100};
-  Double_t xminPadEqual[6] = { 0.0, 0.0, -0.5,    0, -250,   0}; 
-  Double_t xmaxPadEqual[6] = { 2.0, 2.0,  3.5, 13000,  250,   3};
-  TString axisNamePadEqual[6]   = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"};
-  TString axisTitlePadEqual[6]  = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength", "1/pt"};
+  Int_t binsPadEqual[5]    = { 400, 400,    4,   10,   10};
+  Double_t xminPadEqual[5] = { 0.0, 0.0, -0.5,    0, -250}; 
+  Double_t xmaxPadEqual[5] = { 2.0, 2.0,  3.5, 13000, 250};
+  TString axisNamePadEqual[5]   = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength"};
+  TString axisTitlePadEqual[5]  = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength"};
   //
-  fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift, 4:1/pt", 6, binsPadEqual, xminPadEqual, xmaxPadEqual);
-  for (Int_t iaxis=0; iaxis<6;iaxis++){
+  fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift", 5, binsPadEqual, xminPadEqual, xmaxPadEqual);
+  for (Int_t iaxis=0; iaxis<5;iaxis++){
     fHistPadEqual->GetAxis(iaxis)->SetName(axisNamePadEqual[iaxis]);
     fHistPadEqual->GetAxis(iaxis)->SetTitle(axisTitlePadEqual[iaxis]);
   }
@@ -226,6 +227,7 @@ AliTPCcalibGainMult::~AliTPCcalibGainMult(){
   delete fHistdEdxMax;
   delete fHistdEdxTot;
   delete fdEdxTree;
+  if (fBBParam) delete fBBParam;
 }
 
 
@@ -277,7 +279,7 @@ void AliTPCcalibGainMult::Process(AliESDEvent *event) {
   //
   ProcessV0s(event);   // 
   ProcessTOF(event);   //
-  ProcessKinks(event); // not relyable
+  //ProcessKinks(event); // not relyable
   DumpHPT(event);      // 
   UInt_t runNumber = event->GetRunNumber();
   Int_t nContributors = event->GetNumberOfTracks();
@@ -322,7 +324,7 @@ void AliTPCcalibGainMult::Process(AliESDEvent *event) {
     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
     }    
-    if (seed) DumpTrack(track, friendTrack,seed,i); // MI implementation for the identified particles
+    //if (seed) DumpTrack(track, friendTrack,seed,i); // MI implementation for the identified particles
     //
     if (seed) { // seed the container with track parameters and the clusters
       // 
@@ -373,8 +375,8 @@ void AliTPCcalibGainMult::Process(AliESDEvent *event) {
       Float_t meanTot = (1/3.)*(signalArrayTot[0] + signalArrayTot[1] + signalArrayTot[2]); 
       if (meanMax < 1e-5 || meanTot < 1e-5) continue;
       for(Int_t ipad = 0; ipad < 4; ipad ++) {
-       Double_t vecPadEqual[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()};
-       fHistPadEqual->Fill(vecPadEqual);
+       Double_t vecPadEqual[5] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift};
+       if ( TMath::Abs(meanP-kMIPPt)<0.05 ) fHistPadEqual->Fill(vecPadEqual);
       }
       //
       //      if (meanP > 0.4 && meanP < 0.55) {
@@ -462,6 +464,7 @@ Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
   // merging of the component
   //
 
+  const Int_t kMaxEntriesSparse=2000000; // MI- temporary - restrict the THnSparse size
   TIterator* iter = li->MakeIterator();
   AliTPCcalibGainMult* cal = 0;
 
@@ -476,7 +479,10 @@ Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
     if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
     if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector());
     if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
-    if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult());
+    if (cal->GetHistGainMult()) {
+       if (fHistGainMult->GetEntries()<kMaxEntriesSparse) fHistGainMult->Add(cal->GetHistGainMult());
+    } 
+
     if (cal->fHistdEdxMap){
       if (fHistdEdxMap) fHistdEdxMap->Add(cal->fHistdEdxMap);
     }
@@ -523,7 +529,7 @@ Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
     //}
     
   }
-  
+  delete iter;
   return 0;
   
 }
@@ -667,8 +673,8 @@ void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftra
     kp5=(*fBBParam)[4];
   }
   //
-  //AliTPCROC *roc = AliTPCROC::Instance();
-  TDatabasePDG *pdg = TDatabasePDG::Instance();
+  static const AliTPCROC *roc = AliTPCROC::Instance();
+  static const TDatabasePDG *pdg = TDatabasePDG::Instance();
 
   Int_t nclITS   = track->GetNcls(0);
   Int_t ncl   = track->GetTPCncls();
@@ -943,7 +949,7 @@ void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftra
   //
   // Fill histograms
   //
-  if ((isKaon||isProton||isPion||isElectron||isMIP||isMuon&&(!isHighPt)) && dedxDef>0) {
+  if (((isKaon||isProton||isPion||isElectron||isMIP||isMuon)&&(!isHighPt)) && dedxDef>0) {
     //
     Int_t ncont = vertex->GetNContributors();
     for (Int_t ipad=0; ipad<3; ipad++){
@@ -986,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;
   //
   //
@@ -1078,9 +1084,10 @@ void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftra
       "pout="<<pout<<
       "tgl="<<tgl<<
       "track.="<<track<<
-      "trackIn.="<<trackIn<<
-      "trackOut.="<<trackOut<<
-      "tpcOut.="<<tpcOut<<
+      //"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<<
@@ -1153,7 +1160,7 @@ void AliTPCcalibGainMult::ProcessV0s(AliESDEvent * event){
   if (esdFriend->TestSkipBit()) return;
   //
   // 
-  TDatabasePDG *pdg = TDatabasePDG::Instance();  
+  static const TDatabasePDG *pdg = TDatabasePDG::Instance();  
   const Double_t kChi2Cut=5;
   const Double_t kMinR=2;
   const Int_t    kMinNcl=110;
@@ -1775,6 +1782,105 @@ 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;
+  
+  if (!fHistGainSector) return NULL;
+  if (!fHistGainSector->GetAxis(2)) return NULL;
+  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