PWG2/FLOW/macros/makeWeights.C

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