Ref storage set in initialisation
[u/mrichter/AliRoot.git] / PWG2 / runProtonCorrection.C
CommitLineData
b5fc8a3e 1Bool_t runProtonCorrection(Int_t stats = 0, const char* dataset = 0x0) {
2 //macro used to extract the correction maps
3 //using the official correction framework of ALICE
4 //for protons and antiprotons
5 //Author: Panos.Christakoglou@cern.ch
6
7 //________________________________________//
8 //Connect to proof
96f84c25 9 //TProof::Reset("alicecaf.cern.ch");
10 TProof::Open("alicecaf.cern.ch");
b5fc8a3e 11
12 // Enable the STEERBase Package
13 gProof->UploadPackage("STEERBase.par");
14 gProof->EnablePackage("STEERBase");
15 // Enable the ESD Package
16 gProof->UploadPackage("ESD.par");
17 gProof->EnablePackage("ESD");
18 // Enable the AOD Package
19 gProof->UploadPackage("AOD.par");
20 gProof->EnablePackage("AOD");
21 // Enable the Analysis Package
22 gProof->UploadPackage("ANALYSIS.par");
23 gProof->EnablePackage("ANALYSIS");
24 gProof->UploadPackage("ANALYSISalice.par");
25 gProof->EnablePackage("ANALYSISalice");
26 // Enable the CORRFW Package
27 gProof->UploadPackage("CORRFW.par");
28 gProof->EnablePackage("CORRFW");
29
30 gProof->Load("./AliProtonCorrectionTask.cxx+g");
31
32 //________________________________________//
33 //Container definition
34 //Variables of the GRID
35 //For the time being: y-pT
36 //Next step: add Vz
37 //Setting up the container grid...
340705b3 38 //===============//
39 //Global tracking//
40 //===============//
b5fc8a3e 41 const Double_t ymin = -1.0;
42 const Double_t ymax = 1.0;
340705b3 43 const Int_t nbin1 = 20; //bins in y
b5fc8a3e 44 const Double_t ptmin = 0.4;
45 const Double_t ptmax = 3.1;
340705b3 46 const Int_t nbin2 = 27; //bins in pT
47 //===============//
48 // TPC tracking //
49 //===============//
50 /*const Double_t ymin = -0.5;
51 const Double_t ymax = 0.5;
52 const Int_t nbin1 = 10; //bins in y
53 const Double_t ptmin = 0.4;
54 const Double_t ptmax = 0.9;
55 const Int_t nbin2 = 15;*/ //bins in pT
56
b5fc8a3e 57 //Setting up the container grid...
58 UInt_t nstep = 4; //number of selection steps MC
340705b3 59 UInt_t iy = 0;
60 UInt_t ipT = 1;
b5fc8a3e 61 const Int_t nvar = 2; //number of variables on the grid:y-pT
340705b3 62 //arrays for the number of bins in each dimension
b5fc8a3e 63 Int_t iBin[nvar];
64 iBin[0]=nbin1;
65 iBin[1]=nbin2;
66 //arrays for lower bounds :
67 Double_t *binLim1=new Double_t[nbin1+1];
68 Double_t *binLim2=new Double_t[nbin2+1];
69 //values for bin lower bounds
70 for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin + (ymax-ymin) /nbin1*(Double_t)i;
71 for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)ptmin + (ptmax-ptmin)/nbin2*(Double_t)i;
3079041e 72 //CF container for protons
251e4034 73 AliCFContainer* containerProtons = new AliCFContainer("containerProtons","container for protons",nstep,nvar,iBin);
b5fc8a3e 74 //setting the bin limits
3079041e 75 containerProtons->SetBinLimits(iy,binLim1);
76 containerProtons->SetBinLimits(ipT,binLim2);
251e4034 77 //CF container for antiprotons
78 AliCFContainer* containerAntiProtons = new AliCFContainer("containerAntiProtons","container for antiprotons",nstep,nvar,iBin);
3079041e 79 //setting the bin limits
80 containerAntiProtons->SetBinLimits(iy,binLim1);
81 containerAntiProtons->SetBinLimits(ipT,binLim2);
b5fc8a3e 82
83 //________________________________________//
84 // SET TLIST FOR QA HISTOS
85 TList* qaList = new TList();
86 //Cuts
87 const Int_t mintrackrefsTPC = 2;
88 const Int_t mintrackrefsITS = 3;
3079041e 89 const Int_t chargeProtons = 1;
90 const Int_t PDGProtons = 2212;
91 const Int_t chargeAntiProtons = -1;
92 const Int_t PDGAntiProtons = -2212;
340705b3 93
b5fc8a3e 94 const Int_t minclustersTPC = 50;
340705b3 95 const Float_t maxChi2PerTPCCluster = 3.5;
96 const Float_t maxCov11 = 2.0;
97 const Float_t maxCov22 = 2.0;
98 const Float_t maxCov33 = 0.5;
99 const Float_t maxCov44 = 0.5;
100 const Float_t maxCov55 = 2.0;
101 const Float_t maxSigmaToVertexTPC = 2.5;
102
103 const Int_t minclustersITS = 5;
104 const Float_t maxSigmaToVertex = 2.5;
105
b5fc8a3e 106 // Gen-Level kinematic cuts
3079041e 107 AliCFTrackKineCuts *mcKineCutsProtons = new AliCFTrackKineCuts("mcKineCutsProtons",
108 "MC-level kinematic cuts");
109 mcKineCutsProtons->SetPtRange(ptmin,ptmax);
110 mcKineCutsProtons->SetRapidityRange(ymin,ymax);
111 mcKineCutsProtons->SetChargeMC(chargeProtons);
112 mcKineCutsProtons->SetQAOn(qaList);
113
114 AliCFTrackKineCuts *mcKineCutsAntiProtons = new AliCFTrackKineCuts("mcKineCutsAntiProtons",
115 "MC-level kinematic cuts");
116 mcKineCutsAntiProtons->SetPtRange(ptmin,ptmax);
117 mcKineCutsAntiProtons->SetRapidityRange(ymin,ymax);
118 mcKineCutsAntiProtons->SetChargeMC(chargeAntiProtons);
119 mcKineCutsAntiProtons->SetQAOn(qaList);
b5fc8a3e 120
121 //Particle-Level cuts:
3079041e 122 AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts",
123 "MC particle generation cuts");
b5fc8a3e 124 mcGenCuts->SetRequireIsPrimary();
3079041e 125 mcGenCuts->SetRequirePdgCode(PDGProtons);
b5fc8a3e 126 mcGenCuts->SetQAOn(qaList);
127
128 //Acceptance Cuts
3079041e 129 AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts",
130 "MC acceptance cuts");
b5fc8a3e 131 mcAccCuts->SetMinNHitITS(mintrackrefsITS);
132 mcAccCuts->SetMinNHitTPC(mintrackrefsTPC);
133 mcAccCuts->SetQAOn(qaList);
134
135 // Rec-Level kinematic cuts
3079041e 136 AliCFTrackKineCuts *recKineCutsProtons = new AliCFTrackKineCuts("recKineCutsProtons",
137 "rec-level kine cuts");
138 recKineCutsProtons->SetPtRange(ptmin,ptmax);
139 recKineCutsProtons->SetRapidityRange(ymin,ymax);
140 recKineCutsProtons->SetChargeRec(chargeProtons);
141 recKineCutsProtons->SetQAOn(qaList);
142
143 AliCFTrackKineCuts *recKineCutsAntiProtons = new AliCFTrackKineCuts("recKineCutsAntiProtons",
144 "rec-level kine cuts");
145 recKineCutsAntiProtons->SetPtRange(ptmin,ptmax);
146 recKineCutsAntiProtons->SetRapidityRange(ymin,ymax);
147 recKineCutsAntiProtons->SetChargeRec(chargeAntiProtons);
148 recKineCutsAntiProtons->SetQAOn(qaList);
b5fc8a3e 149
3079041e 150 AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts",
151 "rec-level quality cuts");
b5fc8a3e 152 recQualityCuts->SetMinNClusterTPC(minclustersTPC);
340705b3 153 recQualityCuts->SetMaxChi2PerClusterTPC(maxChi2PerTPCCluster);
154 recQualityCuts->SetMaxCovDiagonalElements(maxCov11,maxCov22,maxCov33,maxCov44,maxCov55);
155 recQualityCuts->SetRequireTPCRefit(kTRUE);
156
157 //recQualityCuts->SetRequireITSRefit(kTRUE);
b5fc8a3e 158 recQualityCuts->SetQAOn(qaList);
159
3079041e 160 AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts",
161 "rec-level isPrimary cuts");
340705b3 162 recIsPrimaryCuts->SetMaxNSigmaToVertex(maxSigmaToVertex);
b5fc8a3e 163 recIsPrimaryCuts->SetQAOn(qaList);
164
165 AliCFTrackCutPid* cutPID = new AliCFTrackCutPid("cutPID","ESD_PID");
166 int n_species = AliPID::kSPECIES;
167 Double_t* prior = new Double_t[n_species];
168
169 prior[0] = 0.0244519;
170 prior[1] = 0.0143988;
171 prior[2] = 0.805747 ;
172 prior[3] = 0.0928785;
173 prior[4] = 0.0625243;
174
175 cutPID->SetPriors(prior);
176 cutPID->SetProbabilityCut(0.0);
177 cutPID->SetDetectors("ITS TPC TOF");
340705b3 178 cutPID->SetParticleType(AliPID::kProton , kTRUE);
b5fc8a3e 179 cutPID->SetQAOn(qaList);
180
181 //________________________________________//
3079041e 182 TObjArray* mcListProtons = new TObjArray(0);
183 mcListProtons->AddLast(mcKineCutsProtons);
184 mcListProtons->AddLast(mcGenCuts);
185 TObjArray* mcListAntiProtons = new TObjArray(0);
186 mcListAntiProtons->AddLast(mcKineCutsAntiProtons);
187 mcListAntiProtons->AddLast(mcGenCuts);
b5fc8a3e 188
189 printf("CREATE ACCEPTANCE CUTS\n");
190 TObjArray* accList = new TObjArray(0);
191 accList->AddLast(mcAccCuts);
192
193 printf("CREATE RECONSTRUCTION CUTS\n");
3079041e 194 TObjArray* recListProtons = new TObjArray(0);
195 recListProtons->AddLast(recKineCutsProtons);
196 recListProtons->AddLast(recQualityCuts);
197 recListProtons->AddLast(recIsPrimaryCuts);
198 TObjArray* recListAntiProtons = new TObjArray(0);
199 recListAntiProtons->AddLast(recKineCutsAntiProtons);
200 recListAntiProtons->AddLast(recQualityCuts);
201 recListAntiProtons->AddLast(recIsPrimaryCuts);
b5fc8a3e 202
203 printf("CREATE PID CUTS\n");
204 TObjArray* fPIDCutList = new TObjArray(0);
205 fPIDCutList->AddLast(cutPID);
206
207 //________________________________________//
208 //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
3079041e 209 AliCFManager* manProtons = new AliCFManager();
210 manProtons->SetParticleContainer(containerProtons);
211 manProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListProtons);
212 manProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
213 manProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListProtons);
214 manProtons->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList);
215
216 AliCFManager* manAntiProtons = new AliCFManager();
217 manAntiProtons->SetParticleContainer(containerAntiProtons);
218 manAntiProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListAntiProtons);
219 manAntiProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
220 manAntiProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListAntiProtons);
221 manAntiProtons->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList);
b5fc8a3e 222
223 //________________________________________//
224 //CREATE THE TASK
225 AliProtonCorrectionTask *task = new AliProtonCorrectionTask("AliProtonCorrectionTask");
3079041e 226 task->SetCFManagerProtons(manProtons); //here is set the CF manager
227 task->SetCFManagerAntiProtons(manAntiProtons); //here is set the CF manager
b5fc8a3e 228 task->SetQAList(qaList);
229
230 //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
231 printf("CREATE ANALYSIS MANAGER\n");
232 // Make the analysis manager
233 AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
234
235 AliMCEventHandler* mcHandler = new AliMCEventHandler();
236 AliESDInputHandler* esdHandler = new AliESDInputHandler();
237 mgr->SetMCtruthEventHandler(mcHandler);
238 mgr->SetInputEventHandler(esdHandler);
239
240 //------ input data ------
8a546c82 241 AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
b5fc8a3e 242
243 // ----- output data -----
244 //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
3079041e 245 AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),
246 AliAnalysisManager::kOutputContainer,"corrections.root");
b5fc8a3e 247 // output TH1I for event counting
3079041e 248 AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),
249 AliAnalysisManager::kOutputContainer,"corrections.root");
250 // output Correction Framework Container (for acceptance & efficiency calculations)
251 AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0",
252 AliCFContainer::Class(),
253 AliAnalysisManager::kOutputContainer,"corrections.root");
b5fc8a3e 254 // output Correction Framework Container (for acceptance & efficiency calculations)
3079041e 255 AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("ccontainer1",
256 AliCFContainer::Class(),
257 AliAnalysisManager::kOutputContainer,"corrections.root");
b5fc8a3e 258 // output QA histograms
3079041e 259 AliAnalysisDataContainer *coutput4 = mgr->CreateContainer("clist0",
260 TList::Class(),
261 AliAnalysisManager::kOutputContainer,"corrections.root");
b5fc8a3e 262
263 mgr->AddTask(task);
264 mgr->ConnectInput(task,0,cinput0);
265 mgr->ConnectOutput(task,0,coutput0);
266 mgr->ConnectOutput(task,1,coutput1);
267 mgr->ConnectOutput(task,2,coutput2);
268 mgr->ConnectOutput(task,3,coutput3);
3079041e 269 mgr->ConnectOutput(task,4,coutput4);
b5fc8a3e 270
271 //________________________________________//
272 if (mgr->InitAnalysis()) {
273 if(dataset)
274 mgr->StartAnalysis("proof",dataset,stats);
275 else {
276 // You should get this macro and the txt file from:
277 // http://aliceinfo.cern.ch/Offline/Analysis/CAF/
278 gROOT->LoadMacro("CreateESDChain.C");
279 TChain* chain = 0x0;
280 chain = CreateESDChain("ESD82XX_30K.txt",stats);
281 chain->SetBranchStatus("*Calo*",0);
282
283 mgr->StartAnalysis("proof",chain);
284 }
285 }
286
287 return kTRUE;
288}
289