Removing compilation warnings and failures
[u/mrichter/AliRoot.git] / 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 }