TPC/macros/CalibratedEdxRegion.C

   1 /*
   2   Macro to make an per region dEdx calibration.
   3   Here we assue we have suufiecint ammount of V) PID selected tracks for the calibration
   4
   5   .x $HOME/rootlogon.C
   6   .L $ALICE_ROOT/TPC/macros/CalibratedEdxRegion.C+
   7
   8  */
   9
  10 #include "TTree.h"
  11 #include "TFile.h"
  12 #include "TH1.h"
  13 #include "TH2.h"
  14 #include "TCut.h"
  15 #include "TStatToolkit.h"
  16 #include "TGraphErrors.h"
  17
  18 #include "AliESDtrack.h"
  19 #include "AliESDv0.h"
  20 #include "AliExternalTrackParam.h"
  21 #include "AliTPCdEdxInfo.h"
  22 #include "AliTPCParamSR.h"
  23
  24 //TTree * trees[4]={treeK0, treeLambda,treeALambda,treeGamma};
  25 TTree * trees[4]={0};
  26 TFile * calibrationFile=0;
  27
  28 void InitTrees(const char * v0file= "/hera/alice/miranov/ExpertQA/data/LHC13c/pass1/V0Selected.root"){
  29   //
  30   //
  31   //
  32   //  const char * v0file = "/hera/alice/miranov/ExpertQA/data/LHC13c/pass1/V0Selected.root";
  33   TFile * f = TFile::Open(v0file);
  34   TTree * treeLambda=(TTree*)f->Get("treeLambda");
  35   TTree * treeALambda=(TTree*)f->Get("treeALambda");
  36   TTree * treeK0=(TTree*)f->Get("treeK0");
  37   TTree * treeGamma=(TTree*)f->Get("treeGamma");
  38   calibrationFile=TFile::Open("dEdxCalibration.root","update");
  39   //
  40   // register BetheBloch parameterization for each type of identified particles
  41   //
  42   AliTPCParam param;
  43   param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
  44   treeK0->SetAlias("dEdxExp0","AliTPCParam::BetheBlochAleph(track0.fIp.P()/massPion,1)");
  45   treeK0->SetAlias("dEdxExp1","AliTPCParam::BetheBlochAleph(track1.fIp.P()/massPion,1)");
  46   treeLambda->SetAlias("dEdxExp0","AliTPCParam::BetheBlochAleph(track0.fIp.P()/massProton,1)");
  47   treeLambda->SetAlias("dEdxExp1","AliTPCParam::BetheBlochAleph(track1.fIp.P()/massPion,1)");
  48   treeALambda->SetAlias("dEdxExp0","AliTPCParam::BetheBlochAleph(track0.fIp.P()/massPion,1)");
  49   treeALambda->SetAlias("dEdxExp1","AliTPCParam::BetheBlochAleph(track1.fIp.P()/massProton,1)");
  50   treeGamma->SetAlias("dEdxExp0","AliTPCParam::BetheBlochAleph(track0.fIp.P()/0.005,1)");
  51   treeGamma->SetAlias("dEdxExp1","AliTPCParam::BetheBlochAleph(track1.fIp.P()/0.005,1)");
  52   //
  53   //
  54   //
  55   AliTPCParamSR paramSR;
  56   TString rROC0=TString::Format("%2.2f", 0.5*(paramSR.GetInnerRadiusLow()+paramSR.GetInnerRadiusUp()));
  57   TString rROC1=TString::Format("%2.2f", 0.5*(paramSR.GetPadRowRadii(36,0)+paramSR.GetPadRowRadii(36,paramSR.GetNRowUp1()-1)));
  58   TString rROC2=TString::Format("%2.2f", 0.5*(paramSR.GetPadRowRadii(36,0)+paramSR.GetPadRowRadii(36,paramSR.GetNRowUp()-1)));
  59   //
  60   //
  61   //
  62   Double_t bz=-5;
  63   trees={treeK0, treeLambda,treeALambda,treeGamma};
  64   for (Int_t itree=0; itree<4;itree++){
  65     trees[itree]->SetAlias("bz","-5");
  66     trees[itree]->SetAlias("phi0ROC0",TString::Format("track0.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC0.Data(),bz));
  67     trees[itree]->SetAlias("phi0ROC1",TString::Format("track0.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC1.Data(),bz));
  68     trees[itree]->SetAlias("phi0ROC2",TString::Format("track0.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC2.Data(),bz));
  69     trees[itree]->SetAlias("phi1ROC0",TString::Format("track1.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC0.Data(),bz));
  70     trees[itree]->SetAlias("phi1ROC1",TString::Format("track1.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC1.Data(),bz));
  71     trees[itree]->SetAlias("phi1ROC2",TString::Format("track1.fIp.GetParameterAtRadius(%s+0,%2.2f,7)",rROC2.Data(),bz));
  72     //
  73     //
  74     trees[itree]->SetAlias("side0","(track0.fIp.fP[3]>0&&track0.fIp.fP[1]>0)+2*(track0.fIp.fP[3]<0&&track0.fIp.fP[1]<0)");  // 0 - both sides, 1- A side, 2- C side
  75     trees[itree]->SetAlias("side1","(track1.fIp.fP[3]>0&&track1.fIp.fP[1]>0)+2*(track1.fIp.fP[3]<0&&track1.fIp.fP[1]<0)");
  76   }
  77   /*
  78     consistency check:
  79     // Phi position extrapolated to the IFC OK
  80     treeK0->Draw("track0.GetParameterAtRadius(83.8,-5,7)-track0.fIp.fAlpha-track0.fIp.fP[0]/track0.fIp.fX>>his(100,-0.02,0.02)","abs(track0.GetParameterAtRadius(83.8,-5,7)-track0.fIp.fAlpha)<0.5","",100000);
  81     // Phi angle extrapolated to the IFC
  82     treeK0->Draw("track0.GetParameterAtRadius(83.8,-5,8)-track0.fIp.Phi()>>his(100,-0.02,0.02)","abs(track0.GetParameterAtRadius(83.8,-5,7)-track0.fIp.fAlpha)<0.5","",100000);
  83     //
  84     treeK0->Draw("track0.GetParameterAtRadius(103.8,-5,9)-(track0.fIp.GetParameterAtRadius(103.8,-5,8)-track0.fIp.GetParameterAtRadius(103.8,-5,7))>>his(100,-0.05,0.05)","abs(track0.GetParameterAtRadius(103.8,-5,8)-track0.GetParameterAtRadius(103.8,-5,7))<0.5","",10000);
  85
  86    */
  87 }
  88
  89 void MakeSectorCalibration(){
  90   //
  91   // Get sector phi correction - separate for
  92   //
  93   TCut cutASide0="side0>0";
  94   TCut cutPt0="abs(track0.fIp.fP[4])<3";
  95   TCut cutASide1="side1>0";
  96   TCut cutPt1="abs(track0.fIp.fP[4])<3";
  97   TObjArray * histoArray = new TObjArray(100);
  98   TObjArray * graphArray = new TObjArray(100);
  99   //
 100   //
 101   //
 102   for (Int_t iregion=0; iregion<3; iregion++){
 103     //
 104     TH2F *hisSectorROC=0;
 105     TH2F *hisSectorROCP=0;
 106     TH2F *hisSectorROCM=0;
 107     TGraphErrors * graphs[3]={0};
 108
 109     for (Int_t itree=0; itree<3; itree++){
 110       TString var0=TString::Format("track0.fTPCdEdxInfo.fTPCsignalRegion[%d]/dEdxExp0:18*(phi0ROC%d/pi+(side0-1))>>hisSectorROCPlus%d(36,0,36,60,20,80)",iregion,iregion,iregion);
 111       TString var1=TString::Format("track1.fTPCdEdxInfo.fTPCsignalRegion[%d]/dEdxExp1:18*(phi1ROC%d/pi+(side1-1))>>hisSectorROCMinus%d(36,0,36,60,20,80)",iregion,iregion,iregion);
 112
 113       trees[itree]->Draw(var0.Data(),cutASide0+cutPt0,"colzgoff");
 114       if (hisSectorROCP==0) {
 115         hisSectorROCP=(TH2F*)trees[itree]->GetHistogram()->Clone();
 116       }
 117       if (hisSectorROCP) hisSectorROCP->Add((TH2F*)trees[itree]->GetHistogram());
 118       trees[itree]->Draw(var1.Data(),cutASide1+cutPt1,"colzgoff");
 119       if (hisSectorROCM==0) hisSectorROCM=(TH2F*)trees[itree]->GetHistogram()->Clone();
 120       if (hisSectorROCM) hisSectorROCM->Add((TH2F*)trees[itree]->GetHistogram());
 121     }
 122     hisSectorROC=(TH2F*)hisSectorROCP->Clone();
 123     hisSectorROC->SetName(TString::Format("hisSectorROCBoth%d(36,0,36,60,20,80)",iregion));
 124     hisSectorROC->Add(hisSectorROCM);
 125     //
 126     graphs[0]=TStatToolkit::MakeStat1D(hisSectorROC, 0, 0.85,4,20,1);
 127     graphs[1]=TStatToolkit::MakeStat1D(hisSectorROCP, 0, 0.85,4,24,2);
 128     graphs[2]=TStatToolkit::MakeStat1D(hisSectorROCM, 0, 0.85,4,25,4);
 129     graphs[0]->SetName(TString::Format("graphSectorROCBoth%d(36,0,36,60,20,80)",iregion));
 130     graphs[1]->SetName(TString::Format("graphSectorROCPlus%d(36,0,36,60,20,80)",iregion));
 131     graphs[2]->SetName(TString::Format("graphSectorROCMinus%d(36,0,36,60,20,80)",iregion));
 132     graphs[0]->Draw("alp");
 133     TStatToolkit::MakeStat1D(hisSectorROCP, 0, 0.85,4,24,2)->Draw("lp");
 134     TStatToolkit::MakeStat1D(hisSectorROCM, 0, 0.85,4,25,4)->Draw("lp");
 135     histoArray->AddLast(hisSectorROC);
 136     histoArray->AddLast(hisSectorROCP);
 137     histoArray->AddLast(hisSectorROCM);
 138     for (Int_t i=0; i<3; i++){
 139       graphArray->AddLast(graphs[i]);
 140     }
 141   }
 142   calibrationFile->mkdir("histos");
 143   calibrationFile->cd("histos");
 144   histoArray->Write();
 145   calibrationFile->mkdir("graphs");
 146   calibrationFile->cd("graphs");
 147   graphArray->Write();
 148
 149 }