1 //============================================================================//
2 // Macro makeWeights.C is used to make phi, pt and eta weights. Before using //
3 // this macro you should already have available output file "AnalysisResults" //
4 // with results from various flow analysis methods. When calling this macro //
5 // you must specify the analysis type and flow analysis method whose output //
6 // file will be used to make the weights for the subsequent runs. See bellow //
7 // for available options for analysis type and flow analysis method. Remark: //
8 // for cumulants, GFC and QC, you must also specify the cumulant order. //
9 //============================================================================//
11 //========================= phi-weights ===================================//
12 // Phi-weights are obtained by inverting the azimuthal acceptance profile. //
13 // This procedure fails if there is a phi bin with no entries. In order to //
14 // avoid this one can attempt to rebin original histogram and make binning //
16 Bool_t rebinOriginalPhiHistogram = kTRUE;
17 Int_t nMerge = 4; // indicates how many original bins will be merged into one new bin
19 //========================= pt-weights =========================
20 // You can make pt-weights in three different ways:
21 // 1.) pt-weights are growing linearly as a function of pt for
22 // pt <= ptCutoff. For pt > ptCutoff the pt-weights are
23 // constant and equal to ptMax. To enable this option set
24 // useLinearPtWeights = kTRUE;
25 // 2.) pt-weights are growing quadratically as a function of pt
26 // for pt <= ptCutoff. For pt > ptCutoff the pt-weights are
27 // constant and equal to ptMax. To enable this option set
28 // useQadraticPtWeights = kTRUE;
29 // 3.) pt-weights are simply v(pt) from the specified method's
30 // result for differential flow vs pt. To enable this option
31 // set useLinearPtWeights = kFALSE and useQadraticPtWeights = kFALSE.
32 Double_t ptCutoff = 2; // in GeV
33 Double_t ptWeightsMax = 0.2;
34 Bool_t useLinearPtWeights = kTRUE;
35 Bool_t useQuadraticPtWeights = kFALSE;
36 //========================= eta-weights ========================
37 // Eta-weights are simply v(eta) from the specified method's
38 // result for differential flow vs eta.
39 //==============================================================
41 // Name of the common output file:
42 TString outputFileName = "AnalysisResults.root";
43 //TString outputFileName = "outputCentrality0.root";
45 enum libModes {mLocal,mLocalSource};
47 //void makeWeights(TString type="", TString method="GFC", TString cumulantOrder="4th", Int_t mode=mLocalSource)
48 void makeWeights(TString type="ESD", TString method="SP", TString cumulantOrder="", Int_t mode=mLocal)
50 // 1. type: "ESD", "AOD", "ESDMCkineESD", "ESDMCkineMC", for Monte Carlo and 'on the fly' use simply "", ;
51 // 2. method: MCEP, LYZ1, LYZ2, LYZEP, SP, FQD, GFC or QC;
52 // 3. cumulantOrder: 2nd, 4th, 6th or 8th;
53 // 4. mode: if mode = mLocal -> analyze data on your computer using aliroot
54 // if mode = mLocalSource -> analyze data on your computer using root + source files
56 // Cross-check if the user's settings make sense:
59 // Load needed libraries:
60 LoadLibrariesMW(mode);
62 // Access the common output file:
63 TFile *outputFile = NULL;
64 if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
66 outputFile = TFile::Open(outputFileName.Data(),"READ");
70 cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl;
71 cout<<" directory "<<gSystem->pwd()<<" !!!!"<<endl;
76 // Access the output file of the specified method in the common output file:
77 TString methodFileName = "output";
78 ((methodFileName.Append(method.Data())).Append("analysis")).Append(type.Data())));
79 TDirectoryFile *methodFile = (TDirectoryFile*)outputFile->FindObjectAny(methodFileName.Data());
80 TList *methodList = NULL;
83 TList* listTemp = methodFile->GetListOfKeys();
84 if(listTemp && listTemp->GetEntries() == 1)
86 TString listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast instead)
87 methodFile->GetObject(listName.Data(),methodList);
90 cout<<" WARNING: Accessing TList from TDirectoryFile failed for method "<<method.Data()<<" !!!!"<<endl;
91 cout<<" Did you actually used "<<method.Data()<<" in the analysis?"<<endl;
96 cout<<" WARNING: Couldn't find a TDirectoryFile "<<methodFileName.Data()<<".root !!!!"<<endl;
99 if(!methodList){cout<<" WARNING: methodList is NULL !!!!"<<endl;exit(0);}
101 // Accessing common control and results histograms from which phi, pt and eta weights will be made:
102 AliFlowCommonHist *commonHist = NULL;
103 AliFlowCommonHistResults *commonHistRes = NULL;
104 if(!(method=="GFC"||method=="QC"))
106 commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%s",method.Data()));
107 commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%s",method.Data()));
108 } else if(method=="GFC") // GFC has distinct common hist results for different cumulant orders (but control histos are the same for different orders)
110 commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%s",method.Data()));
111 commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%sOrder%s",cumulantOrder.Data(),method.Data()));
112 } else // this is for sure QC - treated separately because it has distinct both common control and result histograms for different QC orders
114 commonHist = dynamic_cast<AliFlowCommonHist*> methodList->FindObject(Form("AliFlowCommonHist%sOrder%s",cumulantOrder.Data(),method.Data()));
115 commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%sOrder%s",cumulantOrder.Data(),method.Data()));
120 cout<<"WARNING: commonHist is NULL !!!!"<<endl;
127 cout<<"WARNING: commonHistRes is NULL !!!!"<<endl;
132 // Making the output file "weights.root" in which list "weights" will be saved.
133 // List "weights" will hold histograms phiWeights, ptWeights and etaWeights which
134 // will hold the phi, pt and eta weights, respectively:
135 TFile *weightsFile = new TFile("weights.root","RECREATE");
136 TList *weightsList = new TList();
137 gStyle->SetOptStat(0); // remove statistic box from all histograms
140 TH1F *phiWeights = (TH1F*)commonHist->GetHistPhiRP()->Clone("phi_weights"); // to be improved (transferred into TH1D eventually)
141 phiWeights->SetTitle("#phi-weights: correcting for non-uniform acceptance");
142 phiWeights->SetYTitle("w_{#phi}");
143 phiWeights->SetXTitle("#phi");
144 if(rebinOriginalPhiHistogram)
146 phiWeights->Rebin(nMerge);
148 Int_t nBinsPhi = 0; // number of phi bins
149 Int_t counterOfEmptyBinsPhi = 0; // number of empty phi bins
150 Double_t nParticlesInBin = 0.; // number of particles in particular phi bin
151 Double_t nParticlesPerBin = 0.; // average number of particles per phi bin
152 Double_t nParticles = 0.; // number of particles in all phi bins
153 // calculate phi-weights:
154 nBinsPhi = phiWeights->GetNbinsX();
155 nParticles = phiWeights->Integral();
156 if(nBinsPhi) nParticlesPerBin = nParticles/nBinsPhi;
157 for(Int_t b=1;b<=nBinsPhi;b++)
159 Double_t wPhi = 0.; // phi-weight for particular phi bin
160 nParticlesInBin = phiWeights->GetBinContent(b);
163 wPhi = nParticlesPerBin/nParticlesInBin;
166 counterOfEmptyBinsPhi++;
168 phiWeights->SetBinContent(b,wPhi);
169 phiWeights->SetBinError(b,0.);
171 if(!counterOfEmptyBinsPhi)
173 weightsList->Add(phiWeights);
174 cout<<"Phi weights created."<<endl;
177 cout<<"WARNING: Couldn't create phi weights because "<<counterOfEmptyBinsPhi<<" phi bins were empty !!!!"<<endl;
180 // phi-weights for eta subevents:
183 TH1F *phiWeightsSub0 = (TH1F*)commonHist->
184 GetHistPhiSub0()->Clone("phi_weights_sub0");
185 phiWeightsSub0->SetTitle("#phi-weights for subevent 0");
186 phiWeightsSub0->SetYTitle("w_{#phi}");
187 phiWeightsSub0->SetXTitle("#phi");
188 if(rebinOriginalPhiHistogram)
190 phiWeightsSub0->Rebin(nMerge);
193 Int_t nBinsPhiSub0 = 0; // number of phi bins
194 Int_t counterOfEmptyBinsPhiSub0 = 0; // number of empty phi bins
195 Double_t nParticlesInBinSub0 = 0.; // number of particles in this phi bin
196 Double_t nParticlesPerBinSub0 = 0.; // average number of particles/bin
197 Double_t nParticlesSub0 = 0.; // number of particles in all phi bins
199 TH1F *phiWeightsSub1 = (TH1F*)commonHist->
200 GetHistPhiSub1()->Clone("phi_weights_sub1");
201 phiWeightsSub1->SetTitle("#phi-weights for subevent 0");
202 phiWeightsSub1->SetYTitle("w_{#phi}");
203 phiWeightsSub1->SetXTitle("#phi");
204 if(rebinOriginalPhiHistogram)
206 phiWeightsSub1->Rebin(nMerge);
209 Int_t nBinsPhiSub1 = 0; // number of phi bins
210 Int_t counterOfEmptyBinsPhiSub1 = 0; // number of empty phi bins
211 Double_t nParticlesInBinSub1 = 0.; // number of particles in this phi bin
212 Double_t nParticlesPerBinSub1 = 0.; // average number of particles/bin
213 Double_t nParticlesSub1 = 0.; // number of particles in all phi bins
215 // calculate phi-weights for subevent 0:
216 nBinsPhiSub0 = phiWeightsSub0->GetNbinsX();
217 nParticlesSub0 = phiWeightsSub0->Integral();
218 if(nBinsPhiSub0) nParticlesPerBinSub0 = nParticlesSub0/nBinsPhiSub0;
219 for(Int_t b=1;b<=nBinsPhiSub0;b++) {
220 Double_t wPhiSub0 = 0.; // phi-weight for particular phi bin
221 nParticlesInBinSub0 = phiWeightsSub0->GetBinContent(b);
222 if(nParticlesInBinSub0) {
223 wPhiSub0 = nParticlesPerBinSub0/nParticlesInBinSub0;
225 counterOfEmptyBinsPhiSub0++;
227 phiWeightsSub0->SetBinContent(b,wPhiSub0);
228 phiWeightsSub0->SetBinError(b,0.);
230 if(!counterOfEmptyBinsPhiSub0) {
231 weightsList->Add(phiWeightsSub0);
232 cout<<"Phi weights created for subevent 0."<<endl;
234 cout<<"WARNING: Couldn't create phi weights for subevent 0 because "<<counterOfEmptyBinsPhiSub0<<" phi bins were empty !!!!"<<endl;
237 // calculate phi-weights for subevent 1:
238 nBinsPhiSub1 = phiWeightsSub1->GetNbinsX();
239 nParticlesSub1 = phiWeightsSub1->Integral();
240 if(nBinsPhiSub1) nParticlesPerBinSub1 = nParticlesSub1/nBinsPhiSub1;
241 for(Int_t b=1;b<=nBinsPhiSub1;b++) {
242 Double_t wPhiSub1 = 0.; // phi-weight for particular phi bin
243 nParticlesInBinSub1 = phiWeightsSub1->GetBinContent(b);
244 if(nParticlesInBinSub1) {
245 wPhiSub1 = nParticlesPerBinSub1/nParticlesInBinSub1;
247 counterOfEmptyBinsPhiSub1++;
249 phiWeightsSub1->SetBinContent(b,wPhiSub1);
250 phiWeightsSub1->SetBinError(b,0.);
252 if(!counterOfEmptyBinsPhiSub1) {
253 weightsList->Add(phiWeightsSub1);
254 cout<<"Phi weights created for subevent 1."<<endl;
256 cout<<"WARNING: Couldn't create phi weights for subevent 1 because "<<counterOfEmptyBinsPhiSub1<<" phi bins were empty !!!!"<<endl;
262 Double_t ptMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
263 Double_t ptMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
264 Int_t nBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
265 Double_t ptBinWidth = 0.;
266 if(nBinsPt) ptBinWidth = (ptMax-ptMin)/nBinsPt;
267 TH1D *ptWeights = new TH1D("pt_weights","",nBinsPt,ptMin,ptMax);
268 ptWeights->SetXTitle("p_{t} [GeV]");
269 ptWeights->SetYTitle("w_{p_{T}}");
270 if(useLinearPtWeights)
272 ptWeights->SetTitle("Linear p_{T}-weights: optimizing the flow signal");
273 } else if(useQuadraticPtWeights)
275 ptWeights->SetTitle("Quadratic p_{T}-weights: optimizing the flow signal");
278 ptWeights->SetTitle("Differential flow as p_{T}-weights: optimizing the flow signal");
280 // calculate pt weights:
281 for(Int_t b=1;b<=nBinsPt;b++)
283 if(useLinearPtWeights)
285 if(ptMin+b*ptBinWidth < ptCutoff)
287 ptWeights->SetBinContent(b,(ptMin+b*ptBinWidth)*(ptWeightsMax/ptCutoff));
290 ptWeights->SetBinContent(b,ptWeightsMax);
294 weightsList->Add(ptWeights);
295 cout<<"Pt weights (linear) created."<<endl;
297 } else if(useQuadraticPtWeights)
299 if(ptMin+b*ptBinWidth < ptCutoff)
301 ptWeights->SetBinContent(b,pow(ptMin+b*ptBinWidth,2.)*(ptWeightsMax/pow(ptCutoff,2.)));
304 ptWeights->SetBinContent(b,ptWeightsMax);
308 weightsList->Add(ptWeights);
309 cout<<"Pt weights (quadratic) created."<<endl;
311 } else // differential flow result is used as a pt-weight:
313 ptWeights->SetBinContent(b,commonHistRes->GetHistDiffFlowPtPOI()->GetBinContent(b));
316 weightsList->Add(ptWeights);
317 cout<<"Pt weights (from differential flow) created."<<endl;
320 } // end of for(Int_t b=1;b<=nBinsPt;b++)
323 TH1D *etaWeights = commonHistRes->GetHistDiffFlowEtaPOI()->Clone("eta_weights");
324 etaWeights->SetXTitle("#eta [GeV]");
325 etaWeights->SetYTitle("w_{#eta}");
326 etaWeights->SetTitle("Differential flow as #eta-weights: optimizing the flow signal");
329 weightsList->Add(etaWeights);
330 cout<<"Eta weights (from differential flow) created."<<endl;
333 // Save list holding histogram with weights:
334 weightsFile->WriteObject(weightsList,"weights","SingleKey");
336 cout<<"New file \"weights.root\" created to hold those phi, pt and eta weights."<<endl;
341 } // end of void makeWeights(TString type="", TString method="QC", TString cumulantOrder="4th", Int_t mode=mLocal)
343 void CrossCheckSettings()
345 // Check in this method if the settings make sense
347 if(useLinearPtWeights && useQuadraticPtWeights)
350 cout<<"WARNING: You cannot set useLinearPtWeights and useQuadraticPtWeights to kTRUE"<<endl;
351 cout<<" at the same time. Please make up your mind."<<endl;
358 cout<<"WARNING: It doesn't make much sense to have ptCutoff < 0."<<endl;
365 cout<<"WARNING: It doesn't make much sense to have ptWeightsMax < 0."<<endl;
370 } // end of void CrossCheckSettings()
372 void LoadLibrariesMW(const libModes mode) {
374 //--------------------------------------
375 // Load the needed libraries most of them already loaded by aliroot
376 //--------------------------------------
377 //gSystem->Load("libTree");
378 gSystem->Load("libGeom");
379 gSystem->Load("libVMC");
380 gSystem->Load("libXMLIO");
381 gSystem->Load("libPhysics");
383 //----------------------------------------------------------
384 // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
385 //----------------------------------------------------------
387 //--------------------------------------------------------
388 // If you want to use already compiled libraries
389 // in the aliroot distribution
390 //--------------------------------------------------------
392 //==================================================================================
393 //load needed libraries:
394 gSystem->AddIncludePath("-I$ROOTSYS/include");
395 //gSystem->Load("libTree");
398 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
399 gSystem->Load("libANALYSIS");
400 gSystem->Load("libPWGflowBase");
401 //cerr<<"libPWGflowBase loaded ..."<<endl;
405 else if (mode==mLocalSource) {
407 // In root inline compile
410 gROOT->LoadMacro("BaseAliFlowCommonConstants.cxx+");
411 gROOT->LoadMacro("BaseAliFlowLYZConstants.cxx+");
414 gROOT->LoadMacro("BaseAliFlowVector.cxx+");
415 gROOT->LoadMacro("BaseAliFlowTrackSimple.cxx+");
416 gROOT->LoadMacro("BaseAliFlowTrackSimpleCuts.cxx+");
417 gROOT->LoadMacro("BaseAliFlowEventSimple.cxx+");
419 // Output histosgrams
420 gROOT->LoadMacro("BaseAliFlowCommonHist.cxx+");
421 gROOT->LoadMacro("BaseAliFlowCommonHistResults.cxx+");
422 gROOT->LoadMacro("BaseAliFlowLYZHist1.cxx+");
423 gROOT->LoadMacro("BaseAliFlowLYZHist2.cxx+");
425 cout << "finished loading macros!" << endl;
429 } // end of void LoadLibrariesMW(const libModes mode)