PWGCF/FLOW/macros/fqd.C

   1 //==========================================================================//
   2 // Macro to refit the q-distribution in FQD method. By using this macro the //
   3 // fitting parameters for q-distribution can be tuned and correspondingly   //
   4 // the estimate for reference flow harmonic obtained from FQD method can be //
   5 // significantly improved.                                                  //
   6 //==========================================================================//
   7
   8 // Fitting parameters for q-distribution:
   9 Double_t vStart = 0.05; // starting value for v fit
  10 Double_t vMin = 0.0; // lower bound for v fit
  11 Double_t vMax = 0.25; // upper bound for v fit
  12 Double_t sigma2Start = 0.75; // starting value for sigma^2 fit
  13 Double_t sigma2Min = 0.5; // lower bound for sigma^2 fit (according to theorists must be >= 0.5)
  14 Double_t sigma2Max = 2.5; // upper bound for sigma^2 fit
  15 Double_t treshold = 5.; // first and last bin taken for fitting are determined as the first and last bin with more than 'treshold' number of entries
  16 Bool_t finalResultIsFromSigma2Fitted = kTRUE; // final saved result is obtained with sigma^2 fitted (kTRUE) or sigma^2 fixed (kFALSE)
  17 Bool_t printOnTheScreen = kTRUE; // print or not the final results on the screen
  18 Bool_t plotResults = kTRUE; // plot on the screen q-distribution and resulting fitting functions (for non-fixed sigma^2 and fixed sigma^2)
  19
  20 // Name of the common output file:
  21 TString outputFileName = "AnalysisResults.root";
  22 //TString outputFileName = "outputCentrality0.root";
  23
  24 enum libModes {mLocal,mLocalSource};
  25
  26 TFile *commonOutputFile = NULL; // common output file "AnalysisResults.root"
  27 TDirectoryFile *dirFileFQD = NULL; // FQD's TDirectoryFile in "AnalysisResults.root"
  28 TList *outputListFQD = NULL; // output list holding all FQD objects
  29
  30 void fqd(TString analysisType="ESD", Int_t analysisMode=mLocal)
  31 {
  32  // 1. analysisType: type of analysis can be ESD, AOD, MC, MK, ....
  33  //                  (if type="" output files are from simulation 'on the fly')
  34  // 2. analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
  35  //                  if analysisMode = mLocalSource -> analyze data on your computer using root + source files
  36
  37  // Cross-check if the user's settings make sense:
  38  CrossCheckSettings();
  39  // Load needed libraries:
  40  LoadLibrariesFQD(analysisMode);
  41  // Get output list which holds all FQD objects:
  42  GetOutputList(analysisType);
  43  // Redo fit of q-distribution:
  44  RedoFit();
  45  // Plot q-distribution and resulting fitting functions:
  46  if(plotResults) Plot();
  47  // Save the new results in the common output file:
  48  Save();
  49
  50 } // end of void fqd(TString type="ESD", Int_t mode=mLocal)
  51
  52 // ===========================================================================================
  53
  54 void RedoFit()
  55 {
  56  // Redo the fit of q-distribution.
  57
  58  AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
  59  fqd->GetOutputHistograms(outputListFQD);
  60  // Set new fitting parameters:
  61  fqd->SetvStart(vStart);
  62  fqd->SetvMin(vMin);
  63  fqd->SetvMax(vMax);
  64  fqd->SetSigma2Start(sigma2Start);
  65  fqd->SetSigma2Min(sigma2Min);
  66  fqd->SetSigma2Max(sigma2Max);
  67  fqd->SetTreshold(treshold);
  68  fqd->SetFinalResultIsFromSigma2Fitted(finalResultIsFromSigma2Fitted);
  69  fqd->SetPrintOnTheScreen(printOnTheScreen);
  70  // Save new fitting parameters:
  71  fqd->StoreFittingParameters();
  72  // Redo fit:
  73  fqd->Finish(kTRUE);
  74
  75 } // end of void RedoFit(TList *outputList)
  76
  77 // ===========================================================================================
  78
  79 void GetOutputList(TString analysisType)
  80 {
  81  // Get output list which holds all FQD objects.
  82
  83  // Access common output file:
  84  commonOutputFile = AccessOutputFile(outputFileName);
  85
  86  // Access from common output file the TDirectoryFile for FQD method
  87  // and from it the output list holding all objects:
  88  GetListWithHistograms(commonOutputFile,analysisType,"FQD");
  89
  90 } // end of TList* GetOutputList(TString analysisType)
  91
  92 // ===========================================================================================
  93
  94 void Plot()
  95 {
  96  // Plot q-distribution and resulting fitting functions.
  97
  98  gROOT->SetStyle("Plain"); // default color is white instead of gray
  99  gStyle->SetOptStat(0); // remove stat. box from all histos
 100  legend = new TLegend(0.6,0.55,0.85,0.7);
 101  legend->SetFillStyle(0);
 102  // q-distribution:
 103  TH1D *qDistribution = dynamic_cast<TH1D*>(outputListFQD->FindObject("fqDistribution"));
 104  Cosmetics(qDistribution);
 105  legend->AddEntry(qDistribution,"q-distribution","f");
 106  qDistribution->Draw();
 107  // resulting fitting functions:
 108  TF1 *fittingFun[2] = {NULL};
 109  TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
 110  Double_t lineColors[2] = {kBlue,kRed};
 111  for(Int_t s2F=0;s2F<2;s2F++)
 112  {
 113   fittingFun[s2F] = dynamic_cast<TF1*>(outputListFQD->FindObject(Form("fFittingFunction, %s",sigmaFlag[s2F].Data())));
 114   fittingFun[s2F]->SetLineColor(lineColors[s2F]);
 115   fittingFun[s2F]->Draw("SAME");
 116   legend->AddEntry(fittingFun[s2F]->GetName(),sigmaFlag[s2F].Data(),"l");
 117  }
 118  legend->Draw("SAME");
 119
 120 } // end of Plot()
 121
 122 // ===========================================================================================
 123
 124 void Save()
 125 {
 126  // Save the new results of fitting for FQD method in the common output file.
 127
 128  dirFileFQD->Add(outputListFQD,kTRUE);
 129  dirFileFQD->Write(dirFileFQD->GetName(),TObject::kSingleKey+TObject::kOverwrite);
 130  //delete commonOutputFile;
 131
 132 } // end of void Save()
 133
 134 // ===========================================================================================
 135
 136 void Cosmetics(TH1D *hist)
 137 {
 138  // Set cosmetics for the q-distribution.
 139
 140  Int_t firstNonEmptyBin = hist->FindFirstBinAbove(0);
 141  Double_t lowRange = hist->GetBinLowEdge(firstNonEmptyBin);
 142  Int_t lastNonEmptyBin = hist->FindLastBinAbove(0);
 143  Double_t upperRange = hist->GetBinLowEdge(lastNonEmptyBin+10);
 144  hist->GetXaxis()->SetRangeUser(lowRange,upperRange);
 145  hist->SetFillColor(16);
 146  hist->SetTitle("Fitted q-distribution");
 147
 148 } // end of void Cosmetics(TH1D *hist)
 149
 150 // ===========================================================================================
 151
 152 TFile* AccessOutputFile(TString outputFileName)
 153 {
 154  // Access the common output file.
 155
 156  TFile *outputFile = NULL;
 157  if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
 158  {
 159   outputFile = TFile::Open(outputFileName.Data(),"UPDATE");
 160  } else
 161    {
 162     cout<<endl;
 163     cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl;
 164     cout<<"         directory "<<gSystem->pwd()<<" !!!!"<<endl;
 165     cout<<endl;
 166     exit(0);
 167    }
 168
 169  return outputFile;
 170
 171 } // end of TFile* AccessOutputFile(TString outputFileName)
 172
 173 // ===========================================================================================
 174
 175 void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
 176 {
 177  // Access from common output file the TDirectoryFile for FQD method and from it
 178  // the output list holding all objects:
 179
 180  TString fileName = "";
 181  TString listName = "";
 182  // Form a file name:
 183  fileName+="output";
 184  fileName+=methodName.Data();
 185  fileName+="analysis";
 186  fileName+=analysisType.Data();
 187  // Access this file:
 188  dirFileFQD = (TDirectoryFile*)outputFile->FindObjectAny(fileName.Data());
 189  // Form a list name for each method:
 190  if(dirFileFQD)
 191  {
 192   TList* listTemp = dirFileFQD->GetListOfKeys();
 193   if(listTemp && listTemp->GetEntries() == 1)
 194   {
 195    listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast)
 196    dirFileFQD->GetObject(listName.Data(),outputListFQD);
 197   } else
 198     {
 199      cout<<" WARNING: Accessing TList from TDirectoryFile failed miserably for FQD method !!!!"<<endl;
 200      cout<<"          Did you actually use FQD in the analysis?"<<endl;
 201      cout<<endl;
 202      exit(0);
 203     }
 204  } else
 205    {
 206     cout<<" WARNING: Couldn't find a TDirectoryFile "<<fileName.Data()<<" in file "<<outputFileName.Data()<<" !!!!"<<endl;
 207     cout<<"          Perhaps wrong analysis type? Can be \"ESD\",\"AOD\",\"\", etc."<<endl;
 208     cout<<"          Or you didn't use FQD in the analysis which have produced "<<outputFileName.Data()<<"?"<<endl;
 209    }
 210
 211 } // end of void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
 212
 213 //===========================================================================================
 214
 215 void CrossCheckSettings()
 216 {
 217  // Check in this method if the settings make sense.
 218
 219  if(sigma2Min<0.5)
 220  {
 221   cout<<endl;
 222   cout<<"WARNING: According to theorists sigma2Min >= 0.5 !!!!"<<endl;
 223   cout<<endl;
 224   exit(0);
 225  }
 226
 227 } // end of void CrossCheckSettings()
 228
 229 //===========================================================================================
 230
 231 void LoadLibrariesFQD(const libModes mode) {
 232
 233   //--------------------------------------
 234   // Load the needed libraries most of them already loaded by aliroot
 235   //--------------------------------------
 236   gSystem->Load("libTree");
 237   gSystem->Load("libGeom");
 238   gSystem->Load("libVMC");
 239   gSystem->Load("libXMLIO");
 240   gSystem->Load("libPhysics");
 241
 242   //----------------------------------------------------------
 243   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 244   //----------------------------------------------------------
 245   if (mode==mLocal) {
 246     //--------------------------------------------------------
 247     // If you want to use already compiled libraries
 248     // in the aliroot distribution
 249     //--------------------------------------------------------
 250
 251   //==================================================================================
 252   //load needed libraries:
 253   gSystem->AddIncludePath("-I$ROOTSYS/include");
 254   gSystem->Load("libTree");
 255
 256   // for AliRoot
 257   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 258   gSystem->Load("libANALYSIS");
 259   gSystem->Load("libPWGflowBase");
 260   cerr<<"libPWGflowBase loaded ..."<<endl;
 261
 262   }
 263
 264   else if (mode==mLocalSource) {
 265
 266     // In root inline compile
 267
 268     // Constants
 269     gROOT->LoadMacro("BaseAliFlowCommonConstants.cxx+");
 270     gROOT->LoadMacro("BaseAliFlowLYZConstants.cxx+");
 271
 272     // Flow event
 273     gROOT->LoadMacro("BaseAliFlowVector.cxx+");
 274     gROOT->LoadMacro("BaseAliFlowTrackSimple.cxx+");
 275     gROOT->LoadMacro("BaseAliFlowTrackSimpleCuts.cxx+");
 276     gROOT->LoadMacro("BaseAliFlowEventSimple.cxx+");
 277
 278     // Output histosgrams
 279     gROOT->LoadMacro("BaseAliFlowCommonHist.cxx+");
 280     gROOT->LoadMacro("BaseAliFlowCommonHistResults.cxx+");
 281     gROOT->LoadMacro("BaseAliFlowLYZHist1.cxx+");
 282     gROOT->LoadMacro("BaseAliFlowLYZHist2.cxx+");
 283
 284     cout << "finished loading macros!" << endl;
 285
 286   } // end of else if (mode==mLocalSource)
 287
 288 } // end of void LoadLibrariesFQD(const libModes mode)