PWG2/FLOW/macros/makeWeights.C

   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 //=======================================================================//
  11
  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 //==============================================================
  39
  40 enum libModes {mLocal,mLocalSource};
  41
  42 //void makeWeights(TString type="", TString method="GFC", TString cumulantOrder="4th", Int_t mode=mLocalSource)
  43 void makeWeights(TString type="ESD", TString method="", TString cumulantOrder="", Int_t mode=mLocal)
  44 {
  45  // 1. type: "ESD", "AOD", "ESDMCkineESD", "ESDMCkineMC", for Monte Carlo and 'on the fly' use simply "", ;
  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
  50
  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";
  59  TFile *outputFile = NULL;
  60  // Access the common output file:
  61  if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
  62  {
  63   outputFile = TFile::Open(outputFileName.Data(),"READ");
  64  } else
  65    {
  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    }
  72
  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)
  79  {
  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    }
  96
  97  // Accessing common control and results histograms from which phi, pt and eta weights will be made:
  98  AliFlowCommonHist *commonHist = NULL;
  99  AliFlowCommonHistResults *commonHistRes = NULL;
 100  if(!(method=="GFC"||method=="QC"))
 101  {
 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)
 105    {
 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
 109      {
 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()));
 112      }
 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  }
 127
 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++)
 150  {
 151   Double_t wPhi = 0.; // phi-weight for particular phi bin
 152   nParticlesInBin = phiWeights->GetBinContent(b);
 153   if(nParticlesInBin)
 154   {
 155    wPhi = nParticlesPerBin/nParticlesInBin;
 156   } else
 157     {
 158      counterOfEmptyBinsPhi++;
 159     }
 160   phiWeights->SetBinContent(b,wPhi);
 161  }
 162  if(!counterOfEmptyBinsPhi)
 163  {
 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    }
 170
 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
 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)
 250  {
 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)
 263   {
 264    if(ptMin+b*ptBinWidth < ptCutoff)
 265    {
 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)
 277     {
 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++)
 300
 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)
 321
 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()
 350
 351 void LoadLibrariesMW(const libModes mode) {
 352
 353   //--------------------------------------
 354   // Load the needed libraries most of them already loaded by aliroot
 355   //--------------------------------------
 356   //gSystem->Load("libTree");
 357   gSystem->Load("libGeom");
 358   gSystem->Load("libVMC");
 359   gSystem->Load("libXMLIO");
 360   gSystem->Load("libPhysics");
 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");
 374     //gSystem->Load("libTree");
 375
 376     // for AliRoot
 377     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 378     gSystem->Load("libANALYSIS");
 379     gSystem->Load("libPWG2flowCommon");
 380     //cerr<<"libPWG2flowCommon loaded ..."<<endl;
 381
 382   }
 383
 384   else if (mode==mLocalSource) {
 385
 386     // In root inline compile
 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+");
 395     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 396     gROOT->LoadMacro("AliFlowCommon/AliFlowEvent.cxx+");
 397     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 398
 399     // Cuts
 400     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 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
 412 } // end of void LoadLibrariesMW(const libModes mode)
 413
 414
 415
 416
 417