PWG2/FLOW/macros/makeWeights.C

   1 //==========================================================================================
   2 // Before using the macro makeWeights.C you should already have available the output .root
   3 // files from various methods from the previous run over data (without any weights). When
   4 // calling this macro you must specify the analysis type and the method from which output
   5 // file you would like to make the weights for the next run (for the cumulants, GFC and QC,
   6 // you must also specify the order):
   7 //
   8 // 1. type of analysis can be: ESD, AOD, MC, ESDMC0 or ESDMC1;
   9 //
  10 // 2. method can be: MCEP, LYZ1, LYZ2, LYZEP, FQD, GFC or QC;
  11 //
  12 // 3. cumulant order can be: 2nd, 4th, 6th or 8th.
  13 //==========================================================================================
  14
  15
  16 enum libModes {mLocal,mLocalSource};
  17 //mLocal: Analyze data on your computer using aliroot
  18 //mLocalSource: Analyze data on your computer using root + source files
  19
  20
  21 //void makeWeights(TString type="", TString method="GFC", TString cumulantOrder="4th", Int_t mode=mLocalSource)
  22 void makeWeights(TString type="ESD", TString method="GFC", TString cumulantOrder="4th", Int_t mode=mLocal)
  23 {
  24  // load needed libraries:
  25   LoadWeightLibraries(mode);
  26
  27
  28  // open the output file from the first run of the specified method:
  29  TString inputFileName = "output";
  30  TFile* file = NULL;
  31  file = TFile::Open(((((inputFileName.Append(method.Data())).Append("analysis")).Append(type.Data())).Append(".root")).Data(), "READ");
  32
  33  // using pt weights linear or quadratic in pt:
  34  Bool_t useLinearPt    = kTRUE;
  35  Bool_t useQuadraticPt = kFALSE;
  36  if(useLinearPt && useQuadraticPt)
  37  {
  38   cout<<" WARNING: you can use pt weights linear or quadratic in pt, but not at the same time. "<<endl;
  39   exit(0);
  40  }
  41  Bool_t useFunctionOfPt = useLinearPt||useQuadraticPt;
  42
  43  // accessing the results:
  44  TString cobj = "cobj";
  45  TString afc  = "AliFlowCommonHist";
  46  TString afcr = "AliFlowCommonHistResults";
  47
  48  TList *pList = NULL;
  49  AliFlowCommonHist *commonHist = NULL;
  50  AliFlowCommonHistResults *commonHistRes = NULL;
  51
  52  if(file)
  53  {
  54   file->GetObject((cobj.Append(method.Data()).Data()),pList);
  55   if(pList)
  56   {
  57    if(!(method=="GFC"||method=="QC"))
  58    {
  59     commonHist    = dynamic_cast<AliFlowCommonHist*> (pList->FindObject((afc.Append(method.Data())).Data()));
  60     commonHistRes = dynamic_cast<AliFlowCommonHistResults*> (pList->FindObject((afcr.Append(method.Data())).Data()));
  61    }else if(method=="GFC")
  62     {
  63      commonHist    = dynamic_cast<AliFlowCommonHist*> (pList->FindObject((afc.Append(method.Data())).Data()));
  64      commonHistRes = dynamic_cast<AliFlowCommonHistResults*> (pList->FindObject((((afcr.Append(cumulantOrder.Data()).Append("Order")).Append(method.Data())).Data())));
  65     }else
  66      {
  67       commonHist    = dynamic_cast<AliFlowCommonHist*> (pList->FindObject(((afc.Append(cumulantOrder.Data())).Append("Order")).Append(method.Data())));
  68       commonHistRes = dynamic_cast<AliFlowCommonHistResults*> (pList->FindObject((((afcr.Append(cumulantOrder.Data()).Append("Order")).Append(method.Data())).Data())));
  69      }
  70   }//end of if(pList)
  71  }//end of if(file)
  72
  73  //making the output file and creating the TList to hold the histograms with weights:
  74  TFile* outputFile = new TFile("weights.root","RECREATE");
  75  TList* listWeights = new TList();
  76  Int_t nBinsPhi = 0;
  77  Double_t nParticlesInBin = 0;  // number of particles in particular phi bin
  78  Double_t nParticlesPerBin = 0; // average number of particles per phi bin
  79  Double_t nParticles = 0;       // number of particles in all phi bins
  80  Double_t wPhi = 0.;
  81
  82  // common control histos:
  83  if(commonHist)
  84  {
  85   // azimuthal acceptance:
  86   (commonHist->GetHistPhiRP())->SetName("phi_weights");
  87   (commonHist->GetHistPhiRP())->SetTitle("phi_weights: correction for non-uniform acceptance");
  88   (commonHist->GetHistPhiRP())->SetYTitle("weights");
  89   (commonHist->GetHistPhiRP())->SetXTitle("#phi");
  90   nBinsPhi = (commonHist->GetHistPhiRP())->GetNbinsX();
  91   nParticles = (commonHist->GetHistPhiRP())->Integral();
  92   if(nBinsPhi) nParticlesPerBin = nParticles/nBinsPhi;
  93   // loop over phi bins:
  94   for(Int_t b=1;b<nBinsPhi+1;b++)
  95   {
  96    nParticlesInBin = (commonHist->GetHistPhiRP())->GetBinContent(b);
  97    // calculate the phi weight wPhi for each bin:
  98    if(nParticlesInBin) wPhi = nParticlesPerBin/nParticlesInBin;
  99    (commonHist->GetHistPhiRP())->SetBinContent(b,wPhi);
 100   }
 101   listWeights->Add(commonHist->GetHistPhiRP());
 102  }else{cout<<" WARNING: the common control histos from the 1st run were not accessed."<<endl;}
 103
 104  // common results histos:
 105  if(commonHistRes)
 106  {
 107   // diff. flow (pt):
 108   (commonHistRes->GetHistDiffFlowPtPOI())->SetName("pt_weights");
 109   if(!useFunctionOfPt) listWeights->Add(commonHistRes->GetHistDiffFlowPtPOI());
 110   // diff. flow (eta):
 111   (commonHistRes->GetHistDiffFlowEtaPOI())->SetName("eta_weights");
 112   listWeights->Add(commonHistRes->GetHistDiffFlowEtaPOI());
 113  }else{cout<<" WARNING: the common results histos from the 1st run were not accessed."<<endl;}
 114
 115  // pt weights linear and quadratic in pt:
 116  if(useFunctionOfPt)
 117  {
 118   Double_t ptMin = AliFlowCommonConstants::GetPtMin();
 119   Double_t ptMax = AliFlowCommonConstants::GetPtMax();
 120   Int_t nBinsPt  = AliFlowCommonConstants::GetNbinsPt();
 121   Double_t ptCutOff = 2.0; // for pt > ptCutOff use constant weights
 122   if(nBinsPt==0)
 123   {
 124    cout<<" WARNING: number of pt bins is 0. "<<endl;
 125    exit(0);
 126   }
 127   Double_t ptBinWidth = 1.*(ptMax-ptMin)/nBinsPt;
 128
 129   TH1D *ptFunction = new TH1D("ptFunction", "ptFunction", nBinsPt, ptMin, ptMax);
 130   ptFunction->SetName("pt_weights");
 131   ptFunction->SetXTitle("p_{t} [GeV]");
 132   ptFunction->SetYTitle("weights");
 133   Int_t power = 0;
 134   if(useLinearPt)
 135   {
 136    power = 1;
 137    ptFunction->SetTitle("p_{t} weights: #propto p_{t}");
 138   }
 139   if(useQuadraticPt)
 140   {
 141    power = 2;
 142    ptFunction->SetTitle("p_{t} weights: #propto p_{t}^{2}");
 143   }
 144   for(Int_t b=1;b<nBinsPt+1;b++)
 145   {
 146    if(ptMin+b*ptBinWidth < ptCutOff)
 147    {
 148     ptFunction->SetBinContent(b, pow(ptMin+b*ptBinWidth, power));
 149    }else
 150     {
 151      ptFunction->SetBinContent(b, pow(ptCutOff, power));
 152     }
 153   }
 154   listWeights->Add(ptFunction);
 155  }
 156
 157  outputFile->WriteObject(listWeights,"weights","SingleKey");
 158
 159  delete listWeights;
 160  delete outputFile;
 161 }
 162
 163
 164 void LoadWeightLibraries(const libModes mode) {
 165
 166   //--------------------------------------
 167   // Load the needed libraries most of them already loaded by aliroot
 168   //--------------------------------------
 169   gSystem->Load("libTree.so");
 170   gSystem->Load("libGeom.so");
 171   gSystem->Load("libVMC.so");
 172   gSystem->Load("libXMLIO.so");
 173   gSystem->Load("libPhysics.so");
 174
 175   //----------------------------------------------------------
 176   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 177   //----------------------------------------------------------
 178   if (mode==mLocal) {
 179     //--------------------------------------------------------
 180     // If you want to use already compiled libraries
 181     // in the aliroot distribution
 182     //--------------------------------------------------------
 183
 184     //==================================================================================
 185     //load needed libraries:
 186     gSystem->AddIncludePath("-I$ROOTSYS/include");
 187     gSystem->Load("libTree.so");
 188
 189     // for AliRoot
 190     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 191     gSystem->Load("libANALYSIS.so");
 192     gSystem->Load("libPWG2flowCommon.so");
 193     cerr<<"libPWG2flowCommon.so loaded ..."<<endl;
 194
 195   }
 196
 197   else if (mode==mLocalSource) {
 198
 199     // In root inline compile
 200
 201    // Constants
 202     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 203     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 204     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 205
 206     // Flow event
 207     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 208     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 209     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 210
 211     // Cuts
 212     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 213
 214     // Output histosgrams
 215     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 216     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 217     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 218     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 219
 220     cout << "finished loading macros!" << endl;
 221
 222   }
 223
 224 }
 225
 226
 227