]>
Commit | Line | Data |
---|---|---|
04c6b875 | 1 | //=======================================================================// |
2 | // Macro makeWeights.C is used to make phi, pt and eta weights. | |
3 | // Before using the macro makeWeights.C you should already have | |
4 | // available the output .root files from various methods stored | |
5 | // in the common output file "AnalysisResults.root". | |
6 | // When calling this macro you must specify the analysis type | |
7 | // and the method whose output file you would like to use | |
8 | // to make the weights for the subsequent runs (for cumulants, GFC and QC, | |
9 | // you must also specify the cumulant order). | |
10 | //=======================================================================// | |
bc1a7705 | 11 | |
2ed70edf | 12 | //========================= phi-weights ======================== |
13 | // Phi-weights are obtained by inverting and normalizing the | |
14 | // azimuthal acceptance profile. This procedure isn't applicable | |
15 | // if the detector has a gap in azimuthal acceptance (i.e. if | |
16 | // there exists a phi bin with no entries in the histogram for | |
17 | // detector's azimuthal acceptance profile). | |
18 | //========================= pt-weights ========================= | |
19 | // You can make pt-weights in three different ways: | |
20 | // 1.) pt-weights are growing linearly as a function of pt for | |
21 | // pt <= ptCutoff. For pt > ptCutoff the pt-weights are | |
22 | // constant and equal to ptMax. To enable this option set | |
23 | // useLinearPtWeights = kTRUE; | |
24 | // 2.) pt-weights are growing quadratically as a function of pt | |
25 | // for pt <= ptCutoff. For pt > ptCutoff the pt-weights are | |
26 | // constant and equal to ptMax. To enable this option set | |
27 | // useQadraticPtWeights = kTRUE; | |
28 | // 3.) pt-weights are simply v(pt) from the specified method's | |
29 | // result for differential flow vs pt. To enable this option | |
30 | // set useLinearPtWeights = kFALSE and useQadraticPtWeights = kFALSE. | |
31 | Double_t ptCutoff = 2; // in GeV | |
32 | Double_t ptWeightsMax = 0.2; | |
33 | Bool_t useLinearPtWeights = kTRUE; | |
34 | Bool_t useQuadraticPtWeights = kFALSE; | |
35 | //========================= eta-weights ======================== | |
36 | // Eta-weights are simply v(eta) from the specified method's | |
37 | // result for differential flow vs eta. | |
38 | //============================================================== | |
b25cc698 | 39 | |
40 | enum libModes {mLocal,mLocalSource}; | |
b25cc698 | 41 | |
42 | //void makeWeights(TString type="", TString method="GFC", TString cumulantOrder="4th", Int_t mode=mLocalSource) | |
04c6b875 | 43 | void makeWeights(TString type="ESD", TString method="", TString cumulantOrder="", Int_t mode=mLocal) |
2ed70edf | 44 | { |
940a5ed1 | 45 | // 1. type: "ESD", "AOD", "ESDMCkineESD", "ESDMCkineMC", for Monte Carlo and 'on the fly' use simply "", ; |
2ed70edf | 46 | // 2. method: MCEP, LYZ1, LYZ2, LYZEP, SP, FQD, GFC or QC; |
47 | // 3. cumulantOrder: 2nd, 4th, 6th or 8th; | |
48 | // 4. mode: if mode = mLocal -> analyze data on your computer using aliroot | |
49 | // if mode = mLocalSource -> analyze data on your computer using root + source files | |
d88d29dc | 50 | |
2ed70edf | 51 | // Cross-check if the user's settings make sense: |
52 | CrossCheckSettings(); | |
53 | ||
54 | // Load needed libraries: | |
55 | LoadLibrariesMW(mode); | |
56 | ||
57 | // Name of the common output file: | |
58 | TString outputFileName = "AnalysisResults.root"; | |
d88d29dc | 59 | TFile *outputFile = NULL; |
2ed70edf | 60 | // Access the common output file: |
61 | if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists))) | |
d88d29dc | 62 | { |
2ed70edf | 63 | outputFile = TFile::Open(outputFileName.Data(),"READ"); |
64 | } else | |
d88d29dc | 65 | { |
2ed70edf | 66 | cout<<endl; |
67 | cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl; | |
68 | cout<<" directory "<<gSystem->pwd()<<" !!!!"<<endl; | |
69 | cout<<endl; | |
70 | exit(0); | |
71 | } | |
d88d29dc | 72 | |
2ed70edf | 73 | // Access the output file of the specified method in the common output file: |
74 | TString methodFileName = "output"; | |
75 | ((methodFileName.Append(method.Data())).Append("analysis")).Append(type.Data()))); | |
76 | TDirectoryFile *methodFile = (TDirectoryFile*)outputFile->FindObjectAny(methodFileName.Data()); | |
77 | TList *methodList = NULL; | |
78 | if(methodFile) | |
6c10242b | 79 | { |
2ed70edf | 80 | methodFile->GetObject(Form("cobj%s",method.Data()),methodList); |
81 | if(!methodList) | |
82 | { | |
83 | cout<<endl; | |
84 | cout<<"WARNING: Couldn't access the list "<<methodList->GetName()<<" !!!!"<<endl; | |
85 | cout<<endl; | |
86 | exit(0); | |
87 | } | |
88 | } else | |
89 | { | |
90 | cout<<endl; | |
91 | cout<<"WARNING: Couldn't find the file "<<methodFileName.Data()<<".root in "<<endl; | |
92 | cout<<" common output file "<<gSystem->pwd()<<"/"<<outputFileName.Data()<<" !!!!"<<endl; | |
93 | cout<<endl; | |
94 | exit(0); | |
95 | } | |
bc1a7705 | 96 | |
2ed70edf | 97 | // Accessing common control and results histograms from which phi, pt and eta weights will be made: |
bc1a7705 | 98 | AliFlowCommonHist *commonHist = NULL; |
99 | AliFlowCommonHistResults *commonHistRes = NULL; | |
2ed70edf | 100 | if(!(method=="GFC"||method=="QC")) |
bc1a7705 | 101 | { |
2ed70edf | 102 | commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%s",method.Data())); |
103 | commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%s",method.Data())); | |
104 | } else if(method=="GFC") // GFC has distinct common hist results for different cumulant orders (but control histos are the same for different orders) | |
bc1a7705 | 105 | { |
2ed70edf | 106 | commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%s",method.Data())); |
107 | commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%sOrder%s",cumulantOrder.Data(),method.Data())); | |
108 | } else // this is for sure QC - treated separately because it has distinct both common control and result histograms for different QC orders | |
bc1a7705 | 109 | { |
2ed70edf | 110 | commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%sOrder%s",cumulantOrder.Data(),method.Data())); |
111 | commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%sOrder%s",cumulantOrder.Data(),method.Data())); | |
bc1a7705 | 112 | } |
2ed70edf | 113 | if(!commonHist) |
114 | { | |
115 | cout<<endl; | |
116 | cout<<"WARNING: commonHist is NULL !!!!"<<endl; | |
117 | cout<<endl; | |
118 | exit(0); | |
119 | } | |
120 | if(!commonHistRes) | |
121 | { | |
122 | cout<<endl; | |
123 | cout<<"WARNING: commonHistRes is NULL !!!!"<<endl; | |
124 | cout<<endl; | |
125 | exit(0); | |
126 | } | |
bc1a7705 | 127 | |
2ed70edf | 128 | // Making the output file "weights.root" in which list "weights" will be saved. |
129 | // List "weights" will hold histograms phiWeights, ptWeights and etaWeights which | |
130 | // will hold the phi, pt and eta weights, respectively: | |
131 | TFile *weightsFile = new TFile("weights.root","RECREATE"); | |
132 | TList *weightsList = new TList(); | |
133 | gStyle->SetOptStat(0); // remove statistic box from all histograms | |
134 | ||
135 | // phi-weights: | |
136 | TH1F *phiWeights = (TH1F*)commonHist->GetHistPhiRP()->Clone("phi_weights"); // to be improved (transferred into TH1D eventually) | |
137 | phiWeights->SetTitle("#phi-weights: correcting for non-uniform acceptance"); | |
138 | phiWeights->SetYTitle("w_{#phi}"); | |
139 | phiWeights->SetXTitle("#phi"); | |
140 | Int_t nBinsPhi = 0; // number of phi bins | |
141 | Int_t counterOfEmptyBinsPhi = 0; // number of empty phi bins | |
142 | Double_t nParticlesInBin = 0.; // number of particles in particular phi bin | |
143 | Double_t nParticlesPerBin = 0.; // average number of particles per phi bin | |
144 | Double_t nParticles = 0.; // number of particles in all phi bins | |
145 | // calculate phi-weights: | |
146 | nBinsPhi = phiWeights->GetNbinsX(); | |
147 | nParticles = phiWeights->Integral(); | |
148 | if(nBinsPhi) nParticlesPerBin = nParticles/nBinsPhi; | |
149 | for(Int_t b=1;b<=nBinsPhi;b++) | |
bc1a7705 | 150 | { |
2ed70edf | 151 | Double_t wPhi = 0.; // phi-weight for particular phi bin |
152 | nParticlesInBin = phiWeights->GetBinContent(b); | |
153 | if(nParticlesInBin) | |
bc1a7705 | 154 | { |
2ed70edf | 155 | wPhi = nParticlesPerBin/nParticlesInBin; |
156 | } else | |
157 | { | |
158 | counterOfEmptyBinsPhi++; | |
159 | } | |
160 | phiWeights->SetBinContent(b,wPhi); | |
161 | } | |
162 | if(!counterOfEmptyBinsPhi) | |
bc1a7705 | 163 | { |
2ed70edf | 164 | weightsList->Add(phiWeights); |
165 | cout<<"Phi weights created."<<endl; | |
166 | } else | |
167 | { | |
168 | cout<<"WARNING: Couldn't create phi weights because "<<counterOfEmptyBinsPhi<<" phi bins were empty !!!!"<<endl; | |
169 | } | |
bc1a7705 | 170 | |
04c6b875 | 171 | // phi-weights for eta subevents: |
172 | if (method=="SP"){ | |
173 | //subevent 0 | |
174 | TH1F *phiWeightsSub0 = (TH1F*)commonHist-> | |
175 | GetHistPhiSub0()->Clone("phi_weights_sub0"); | |
176 | phiWeightsSub0->SetTitle("#phi-weights for subevent 0"); | |
177 | phiWeightsSub0->SetYTitle("w_{#phi}"); | |
178 | phiWeightsSub0->SetXTitle("#phi"); | |
179 | Int_t nBinsPhiSub0 = 0; // number of phi bins | |
180 | Int_t counterOfEmptyBinsPhiSub0 = 0; // number of empty phi bins | |
181 | Double_t nParticlesInBinSub0 = 0.; // number of particles in this phi bin | |
182 | Double_t nParticlesPerBinSub0 = 0.; // average number of particles/bin | |
183 | Double_t nParticlesSub0 = 0.; // number of particles in all phi bins | |
184 | //subevent 1 | |
185 | TH1F *phiWeightsSub1 = (TH1F*)commonHist-> | |
186 | GetHistPhiSub1()->Clone("phi_weights_sub1"); | |
187 | phiWeightsSub1->SetTitle("#phi-weights for subevent 0"); | |
188 | phiWeightsSub1->SetYTitle("w_{#phi}"); | |
189 | phiWeightsSub1->SetXTitle("#phi"); | |
190 | Int_t nBinsPhiSub1 = 0; // number of phi bins | |
191 | Int_t counterOfEmptyBinsPhiSub1 = 0; // number of empty phi bins | |
192 | Double_t nParticlesInBinSub1 = 0.; // number of particles in this phi bin | |
193 | Double_t nParticlesPerBinSub1 = 0.; // average number of particles/bin | |
194 | Double_t nParticlesSub1 = 0.; // number of particles in all phi bins | |
195 | ||
196 | // calculate phi-weights for subevent 0: | |
197 | nBinsPhiSub0 = phiWeightsSub0->GetNbinsX(); | |
198 | nParticlesSub0 = phiWeightsSub0->Integral(); | |
199 | if(nBinsPhiSub0) nParticlesPerBinSub0 = nParticlesSub0/nBinsPhiSub0; | |
200 | for(Int_t b=1;b<=nBinsPhiSub0;b++) { | |
201 | Double_t wPhiSub0 = 0.; // phi-weight for particular phi bin | |
202 | nParticlesInBinSub0 = phiWeightsSub0->GetBinContent(b); | |
203 | if(nParticlesInBinSub0) { | |
204 | wPhiSub0 = nParticlesPerBinSub0/nParticlesInBinSub0; | |
205 | } else { | |
206 | counterOfEmptyBinsPhiSub0++; | |
207 | } | |
208 | phiWeightsSub0->SetBinContent(b,wPhiSub0); | |
209 | } | |
210 | if(!counterOfEmptyBinsPhiSub0) { | |
211 | weightsList->Add(phiWeightsSub0); | |
212 | cout<<"Phi weights created."<<endl; | |
213 | } else { | |
214 | cout<<"WARNING: Couldn't create phi weights for subevent 0 because "<<counterOfEmptyBinsPhiSub0<<" phi bins were empty !!!!"<<endl; | |
215 | } | |
216 | ||
217 | // calculate phi-weights for subevent 1: | |
218 | nBinsPhiSub1 = phiWeightsSub1->GetNbinsX(); | |
219 | nParticlesSub1 = phiWeightsSub1->Integral(); | |
220 | if(nBinsPhiSub1) nParticlesPerBinSub1 = nParticlesSub1/nBinsPhiSub1; | |
221 | for(Int_t b=1;b<=nBinsPhiSub1;b++) { | |
222 | Double_t wPhiSub1 = 0.; // phi-weight for particular phi bin | |
223 | nParticlesInBinSub1 = phiWeightsSub1->GetBinContent(b); | |
224 | if(nParticlesInBinSub1) { | |
225 | wPhiSub1 = nParticlesPerBinSub1/nParticlesInBinSub1; | |
226 | } else { | |
227 | counterOfEmptyBinsPhiSub1++; | |
228 | } | |
229 | phiWeightsSub1->SetBinContent(b,wPhiSub1); | |
230 | } | |
231 | if(!counterOfEmptyBinsPhiSub1) { | |
232 | weightsList->Add(phiWeightsSub1); | |
233 | cout<<"Phi weights created."<<endl; | |
234 | } else { | |
235 | cout<<"WARNING: Couldn't create phi weights for subevent 1 because "<<counterOfEmptyBinsPhiSub1<<" phi bins were empty !!!!"<<endl; | |
236 | } | |
237 | ||
238 | } | |
239 | ||
2ed70edf | 240 | // pt-weights: |
241 | Double_t ptMin = AliFlowCommonConstants::GetMaster()->GetPtMin(); | |
242 | Double_t ptMax = AliFlowCommonConstants::GetMaster()->GetPtMax(); | |
243 | Int_t nBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt(); | |
244 | Double_t ptBinWidth = 0.; | |
245 | if(nBinsPt) ptBinWidth = (ptMax-ptMin)/nBinsPt; | |
246 | TH1D *ptWeights = new TH1D("pt_weights","",nBinsPt,ptMin,ptMax); | |
247 | ptWeights->SetXTitle("p_{t} [GeV]"); | |
248 | ptWeights->SetYTitle("w_{p_{T}}"); | |
249 | if(useLinearPtWeights) | |
6c10242b | 250 | { |
2ed70edf | 251 | ptWeights->SetTitle("Linear p_{T}-weights: optimizing the flow signal"); |
252 | } else if(useQuadraticPtWeights) | |
253 | { | |
254 | ptWeights->SetTitle("Quadratic p_{T}-weights: optimizing the flow signal"); | |
255 | } else | |
256 | { | |
257 | ptWeights->SetTitle("Differential flow as p_{T}-weights: optimizing the flow signal"); | |
258 | } | |
259 | // calculate pt weights: | |
260 | for(Int_t b=1;b<=nBinsPt;b++) | |
261 | { | |
262 | if(useLinearPtWeights) | |
6c10242b | 263 | { |
2ed70edf | 264 | if(ptMin+b*ptBinWidth < ptCutoff) |
6c10242b | 265 | { |
2ed70edf | 266 | ptWeights->SetBinContent(b,(ptMin+b*ptBinWidth)*(ptWeightsMax/ptCutoff)); |
267 | } else | |
268 | { | |
269 | ptWeights->SetBinContent(b,ptWeightsMax); | |
270 | } | |
271 | if(b==nBinsPt) | |
272 | { | |
273 | weightsList->Add(ptWeights); | |
274 | cout<<"Pt weights (linear) created."<<endl; | |
275 | } | |
276 | } else if(useQuadraticPtWeights) | |
6c10242b | 277 | { |
2ed70edf | 278 | if(ptMin+b*ptBinWidth < ptCutoff) |
279 | { | |
280 | ptWeights->SetBinContent(b,pow(ptMin+b*ptBinWidth,2.)*(ptWeightsMax/pow(ptCutoff,2.))); | |
281 | } else | |
282 | { | |
283 | ptWeights->SetBinContent(b,ptWeightsMax); | |
284 | } | |
285 | if(b==nBinsPt) | |
286 | { | |
287 | weightsList->Add(ptWeights); | |
288 | cout<<"Pt weights (quadratic) created."<<endl; | |
289 | } | |
290 | } else // differential flow result is used as a pt-weight: | |
291 | { | |
292 | ptWeights->SetBinContent(b,commonHistRes->GetHistDiffFlowPtPOI()->GetBinContent(b)); | |
293 | if(b==nBinsPt) | |
294 | { | |
295 | weightsList->Add(ptWeights); | |
296 | cout<<"Pt weights (from differential flow) created."<<endl; | |
297 | } | |
298 | } | |
299 | } // end of for(Int_t b=1;b<=nBinsPt;b++) | |
03a02aca | 300 | |
2ed70edf | 301 | // eta-weights: |
302 | TH1D *etaWeights = commonHistRes->GetHistDiffFlowEtaPOI()->Clone("eta_weights"); | |
303 | etaWeights->SetXTitle("#eta [GeV]"); | |
304 | etaWeights->SetYTitle("w_{#eta}"); | |
305 | etaWeights->SetTitle("Differential flow as #eta-weights: optimizing the flow signal"); | |
306 | if(etaWeights) | |
307 | { | |
308 | weightsList->Add(etaWeights); | |
309 | cout<<"Eta weights (from differential flow) created."<<endl; | |
310 | } | |
311 | ||
312 | // Save list holding histogram with weights: | |
313 | weightsFile->WriteObject(weightsList,"weights","SingleKey"); | |
314 | ||
315 | cout<<"New file \"weights.root\" created to hold those phi, pt and eta weights."<<endl; | |
316 | ||
317 | delete weightsList; | |
318 | delete weightsFile; | |
319 | ||
320 | } // end of void makeWeights(TString type="", TString method="QC", TString cumulantOrder="4th", Int_t mode=mLocal) | |
03a02aca | 321 | |
2ed70edf | 322 | void CrossCheckSettings() |
323 | { | |
324 | // Check in this method if the settings make sense | |
325 | ||
326 | if(useLinearPtWeights && useQuadraticPtWeights) | |
327 | { | |
328 | cout<<endl; | |
329 | cout<<"WARNING: You cannot set useLinearPtWeights and useQuadraticPtWeights to kTRUE"<<endl; | |
330 | cout<<" at the same time. Please make up your mind."<<endl; | |
331 | cout<<endl; | |
332 | exit(0); | |
333 | } | |
334 | if(ptCutoff<0.) | |
335 | { | |
336 | cout<<endl; | |
337 | cout<<"WARNING: It doesn't make much sense to have ptCutoff < 0."<<endl; | |
338 | cout<<endl; | |
339 | exit(0); | |
340 | } | |
341 | if(ptWeightsMax<0.) | |
342 | { | |
343 | cout<<endl; | |
344 | cout<<"WARNING: It doesn't make much sense to have ptWeightsMax < 0."<<endl; | |
345 | cout<<endl; | |
346 | exit(0); | |
347 | } | |
348 | ||
349 | } // end of void CrossCheckSettings() | |
03a02aca | 350 | |
2ed70edf | 351 | void LoadLibrariesMW(const libModes mode) { |
b25cc698 | 352 | |
353 | //-------------------------------------- | |
354 | // Load the needed libraries most of them already loaded by aliroot | |
355 | //-------------------------------------- | |
d88d29dc | 356 | //gSystem->Load("libTree"); |
5d040cf3 | 357 | gSystem->Load("libGeom"); |
358 | gSystem->Load("libVMC"); | |
359 | gSystem->Load("libXMLIO"); | |
360 | gSystem->Load("libPhysics"); | |
b25cc698 | 361 | |
362 | //---------------------------------------------------------- | |
363 | // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< | |
364 | //---------------------------------------------------------- | |
365 | if (mode==mLocal) { | |
366 | //-------------------------------------------------------- | |
367 | // If you want to use already compiled libraries | |
368 | // in the aliroot distribution | |
369 | //-------------------------------------------------------- | |
370 | ||
371 | //================================================================================== | |
372 | //load needed libraries: | |
373 | gSystem->AddIncludePath("-I$ROOTSYS/include"); | |
d88d29dc | 374 | //gSystem->Load("libTree"); |
b25cc698 | 375 | |
376 | // for AliRoot | |
377 | gSystem->AddIncludePath("-I$ALICE_ROOT/include"); | |
5d040cf3 | 378 | gSystem->Load("libANALYSIS"); |
379 | gSystem->Load("libPWG2flowCommon"); | |
d88d29dc | 380 | //cerr<<"libPWG2flowCommon loaded ..."<<endl; |
b25cc698 | 381 | |
382 | } | |
383 | ||
384 | else if (mode==mLocalSource) { | |
385 | ||
386 | // In root inline compile | |
d90df6c3 | 387 | |
388 | // Constants | |
389 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+"); | |
390 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+"); | |
391 | gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+"); | |
392 | ||
393 | // Flow event | |
394 | gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); | |
7382279b | 395 | gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+"); |
396 | gROOT->LoadMacro("AliFlowCommon/AliFlowEvent.cxx+"); | |
d90df6c3 | 397 | gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+"); |
398 | ||
399 | // Cuts | |
400 | gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+"); | |
b25cc698 | 401 | |
402 | // Output histosgrams | |
403 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+"); | |
404 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+"); | |
405 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+"); | |
406 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+"); | |
407 | ||
408 | cout << "finished loading macros!" << endl; | |
409 | ||
410 | } | |
411 | ||
2ed70edf | 412 | } // end of void LoadLibrariesMW(const libModes mode) |
413 | ||
414 | ||
b25cc698 | 415 | |
416 | ||
03a02aca | 417 |