PWGCF/FLOW/macros/makeWeights.C

   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 //============================================================================//
  10
  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 //
  15 // in phi coarser:                                                         //
  16 Bool_t rebinOriginalPhiHistogram = kTRUE;
  17 Int_t nMerge = 4; // indicates how many original bins will be merged into one new bin
  18
  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 //==============================================================
  40
  41 // Name of the common output file:
  42 TString outputFileName = "AnalysisResults.root";
  43 //TString outputFileName = "outputCentrality0.root";
  44
  45 enum libModes {mLocal,mLocalSource};
  46
  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)
  49 {
  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
  55
  56  // Cross-check if the user's settings make sense:
  57  CrossCheckSettings();
  58
  59  // Load needed libraries:
  60  LoadLibrariesMW(mode);
  61
  62  // Access the common output file:
  63  TFile *outputFile = NULL;
  64  if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
  65  {
  66   outputFile = TFile::Open(outputFileName.Data(),"READ");
  67  } else
  68    {
  69     cout<<endl;
  70     cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl;
  71     cout<<"         directory "<<gSystem->pwd()<<" !!!!"<<endl;
  72     cout<<endl;
  73     exit(0);
  74    }
  75
  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;
  81  if(methodFile)
  82  {
  83   TList* listTemp = methodFile->GetListOfKeys();
  84   if(listTemp && listTemp->GetEntries() == 1)
  85   {
  86    TString listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast instead)
  87    methodFile->GetObject(listName.Data(),methodList);
  88   } else
  89      {
  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;
  92       cout<<endl;
  93      }
  94  } else
  95    {
  96     cout<<" WARNING: Couldn't find a TDirectoryFile "<<methodFileName.Data()<<".root !!!!"<<endl;
  97    }
  98
  99  if(!methodList){cout<<" WARNING: methodList is NULL !!!!"<<endl;exit(0);}
 100
 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"))
 105  {
 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)
 109    {
 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
 113      {
 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()));
 116      }
 117  if(!commonHist)
 118  {
 119   cout<<endl;
 120   cout<<"WARNING: commonHist is NULL !!!!"<<endl;
 121   cout<<endl;
 122   exit(0);
 123  }
 124  if(!commonHistRes)
 125  {
 126   cout<<endl;
 127   cout<<"WARNING: commonHistRes is NULL !!!!"<<endl;
 128   cout<<endl;
 129   exit(0);
 130  }
 131
 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
 138
 139  // phi-weights:
 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)
 145    {
 146      phiWeights->Rebin(nMerge);
 147    }
 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++)
 158  {
 159   Double_t wPhi = 0.; // phi-weight for particular phi bin
 160   nParticlesInBin = phiWeights->GetBinContent(b);
 161   if(nParticlesInBin)
 162   {
 163    wPhi = nParticlesPerBin/nParticlesInBin;
 164   } else
 165     {
 166      counterOfEmptyBinsPhi++;
 167     }
 168   phiWeights->SetBinContent(b,wPhi);
 169   phiWeights->SetBinError(b,0.);
 170  }
 171  if(!counterOfEmptyBinsPhi)
 172  {
 173   weightsList->Add(phiWeights);
 174   cout<<"Phi weights created."<<endl;
 175  } else
 176    {
 177     cout<<"WARNING: Couldn't create phi weights because "<<counterOfEmptyBinsPhi<<" phi bins were empty !!!!"<<endl;
 178    }
 179
 180 // phi-weights for eta subevents:
 181 if (method=="SP"){
 182   //subevent 0
 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)
 189     {
 190       phiWeightsSub0->Rebin(nMerge);
 191     }
 192
 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
 198   //subevent 1
 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)
 205     {
 206       phiWeightsSub1->Rebin(nMerge);
 207     }
 208
 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
 214
 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;
 224     } else {
 225       counterOfEmptyBinsPhiSub0++;
 226     }
 227     phiWeightsSub0->SetBinContent(b,wPhiSub0);
 228     phiWeightsSub0->SetBinError(b,0.);
 229   }
 230   if(!counterOfEmptyBinsPhiSub0) {
 231     weightsList->Add(phiWeightsSub0);
 232     cout<<"Phi weights created for subevent 0."<<endl;
 233   } else {
 234     cout<<"WARNING: Couldn't create phi weights for subevent 0 because "<<counterOfEmptyBinsPhiSub0<<" phi bins were empty !!!!"<<endl;
 235   }
 236
 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;
 246     } else {
 247       counterOfEmptyBinsPhiSub1++;
 248     }
 249     phiWeightsSub1->SetBinContent(b,wPhiSub1);
 250     phiWeightsSub1->SetBinError(b,0.);
 251   }
 252   if(!counterOfEmptyBinsPhiSub1) {
 253     weightsList->Add(phiWeightsSub1);
 254     cout<<"Phi weights created for subevent 1."<<endl;
 255   } else {
 256     cout<<"WARNING: Couldn't create phi weights for subevent 1 because "<<counterOfEmptyBinsPhiSub1<<" phi bins were empty !!!!"<<endl;
 257   }
 258
 259  }
 260
 261  // pt-weights:
 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)
 271  {
 272   ptWeights->SetTitle("Linear p_{T}-weights: optimizing the flow signal");
 273  } else if(useQuadraticPtWeights)
 274    {
 275     ptWeights->SetTitle("Quadratic p_{T}-weights: optimizing the flow signal");
 276    } else
 277      {
 278       ptWeights->SetTitle("Differential flow as p_{T}-weights: optimizing the flow signal");
 279      }
 280  // calculate pt weights:
 281  for(Int_t b=1;b<=nBinsPt;b++)
 282  {
 283   if(useLinearPtWeights)
 284   {
 285    if(ptMin+b*ptBinWidth < ptCutoff)
 286    {
 287     ptWeights->SetBinContent(b,(ptMin+b*ptBinWidth)*(ptWeightsMax/ptCutoff));
 288    } else
 289      {
 290       ptWeights->SetBinContent(b,ptWeightsMax);
 291      }
 292      if(b==nBinsPt)
 293      {
 294       weightsList->Add(ptWeights);
 295       cout<<"Pt weights (linear) created."<<endl;
 296      }
 297   } else if(useQuadraticPtWeights)
 298     {
 299      if(ptMin+b*ptBinWidth < ptCutoff)
 300      {
 301       ptWeights->SetBinContent(b,pow(ptMin+b*ptBinWidth,2.)*(ptWeightsMax/pow(ptCutoff,2.)));
 302      } else
 303        {
 304         ptWeights->SetBinContent(b,ptWeightsMax);
 305        }
 306        if(b==nBinsPt)
 307        {
 308         weightsList->Add(ptWeights);
 309         cout<<"Pt weights (quadratic) created."<<endl;
 310        }
 311     } else // differential flow result is used as a pt-weight:
 312       {
 313        ptWeights->SetBinContent(b,commonHistRes->GetHistDiffFlowPtPOI()->GetBinContent(b));
 314        if(b==nBinsPt)
 315        {
 316         weightsList->Add(ptWeights);
 317         cout<<"Pt weights (from differential flow) created."<<endl;
 318        }
 319       }
 320  } // end of for(Int_t b=1;b<=nBinsPt;b++)
 321
 322  // eta-weights:
 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");
 327  if(etaWeights)
 328  {
 329   weightsList->Add(etaWeights);
 330   cout<<"Eta weights (from differential flow) created."<<endl;
 331  }
 332
 333  // Save list holding histogram with weights:
 334  weightsFile->WriteObject(weightsList,"weights","SingleKey");
 335
 336  cout<<"New file \"weights.root\" created to hold those phi, pt and eta weights."<<endl;
 337
 338  delete weightsList;
 339  delete weightsFile;
 340
 341 } // end of void makeWeights(TString type="", TString method="QC", TString cumulantOrder="4th", Int_t mode=mLocal)
 342
 343 void CrossCheckSettings()
 344 {
 345  // Check in this method if the settings make sense
 346
 347  if(useLinearPtWeights && useQuadraticPtWeights)
 348  {
 349   cout<<endl;
 350   cout<<"WARNING: You cannot set useLinearPtWeights and useQuadraticPtWeights to kTRUE"<<endl;
 351   cout<<"          at the same time. Please make up your mind."<<endl;
 352   cout<<endl;
 353   exit(0);
 354  }
 355  if(ptCutoff<0.)
 356  {
 357   cout<<endl;
 358   cout<<"WARNING: It doesn't make much sense to have ptCutoff < 0."<<endl;
 359   cout<<endl;
 360   exit(0);
 361  }
 362  if(ptWeightsMax<0.)
 363  {
 364   cout<<endl;
 365   cout<<"WARNING: It doesn't make much sense to have ptWeightsMax < 0."<<endl;
 366   cout<<endl;
 367   exit(0);
 368  }
 369
 370 } // end of void CrossCheckSettings()
 371
 372 void LoadLibrariesMW(const libModes mode) {
 373
 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");
 382
 383   //----------------------------------------------------------
 384   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 385   //----------------------------------------------------------
 386   if (mode==mLocal) {
 387     //--------------------------------------------------------
 388     // If you want to use already compiled libraries
 389     // in the aliroot distribution
 390     //--------------------------------------------------------
 391
 392     //==================================================================================
 393     //load needed libraries:
 394     gSystem->AddIncludePath("-I$ROOTSYS/include");
 395     //gSystem->Load("libTree");
 396
 397     // for AliRoot
 398     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 399     gSystem->Load("libANALYSIS");
 400     gSystem->Load("libPWGflowBase");
 401     //cerr<<"libPWGflowBase loaded ..."<<endl;
 402
 403   }
 404
 405   else if (mode==mLocalSource) {
 406
 407     // In root inline compile
 408
 409    // Constants
 410     gROOT->LoadMacro("BaseAliFlowCommonConstants.cxx+");
 411     gROOT->LoadMacro("BaseAliFlowLYZConstants.cxx+");
 412
 413     // Flow event
 414     gROOT->LoadMacro("BaseAliFlowVector.cxx+");
 415     gROOT->LoadMacro("BaseAliFlowTrackSimple.cxx+");
 416     gROOT->LoadMacro("BaseAliFlowTrackSimpleCuts.cxx+");
 417     gROOT->LoadMacro("BaseAliFlowEventSimple.cxx+");
 418
 419     // Output histosgrams
 420     gROOT->LoadMacro("BaseAliFlowCommonHist.cxx+");
 421     gROOT->LoadMacro("BaseAliFlowCommonHistResults.cxx+");
 422     gROOT->LoadMacro("BaseAliFlowLYZHist1.cxx+");
 423     gROOT->LoadMacro("BaseAliFlowLYZHist2.cxx+");
 424
 425     cout << "finished loading macros!" << endl;
 426
 427   }
 428
 429 } // end of void LoadLibrariesMW(const libModes mode)
 430
 431
 432
 433
 434