PWGCF/FLOW/macros/showSpread.C

   1 // ...
   2
   3 enum libModes {mLocal,mLocalSource};
   4
   5 const Int_t nFilesMax = 10; // number of files to be accessed to estimate spread
   6
   7 Bool_t showPlotForReferenceFlow = kTRUE;
   8 Bool_t showPlotForIntegratedFlowRP = kTRUE;
   9 Bool_t showPlotForIntegratedFlowPOI = kTRUE;
  10 Bool_t showErrorOnMergedResult = kTRUE;
  11
  12 const Int_t nMethods = 13;
  13 TString method[nMethods] = {"MCEP","SP","2,GFC","2,QC","4,GFC","4,QC","6,GFC","6,QC","8,GFC","8,QC","FQD","LYZ1SUM","LYZ1PROD"};
  14 Int_t methodMarkerStyle[nMethods] = {21,21,21,21,21,21,21,21,21,21,21,21,21};
  15 Int_t methodMarkerColor[nMethods] = {kGray+1,kViolet-8,kBlue-9,kRed-7,kBlue-9,kRed-7,kBlue-9,kRed-7,kBlue-9,kRed-7,kOrange-8,kYellow-5,kYellow-2};
  16 Int_t methodMeshColor[nMethods] = {kGray,kViolet-9,kBlue-10,kRed-10,kBlue-10,kRed-10,kBlue-10,kRed-10,kBlue-10,kRed-10,kOrange-9,kYellow-8,kYellow-6};
  17
  18 /*
  19 const Int_t nMethods = 4;
  20 TString method[nMethods] = {"2,QC","4,QC","6,QC","8,QC"};
  21 Int_t methodMarkerStyle[nMethods] = {21,21,21,21};
  22 Int_t methodMarkerColor[nMethods] = {kRed-7,kRed-7,kRed-7,kRed-7};
  23 Int_t methodMeshColor[nMethods] = {kRed-10,kRed-10,kRed-10,kRed-10};
  24 */
  25
  26 void showSpread(TString type="", Int_t mode=mLocal)
  27 {
  28  // type: type of analysis can be ESD, AOD, MC, ESDMCkineESD, ESDMCkineMC
  29  //       (if type="" output files are from MC simulation (default))
  30  // mode: if mode = mLocal: analyze data on your computer using aliroot
  31  //       if mode = mLocalSource: analyze data on your computer using root + source files
  32
  33  // Cross-check if the user's settings make sense:
  34  CrossCheckUserSettings();
  35
  36  // Load needed libraries:
  37  LoadLibrariesSS(mode);
  38
  39  // Output file name:
  40  TString outputFileName = "AnalysisResults.root";
  41
  42  // Labels for reference flow, integrated flow of RPs and of POIs:
  43  TString label[3] = {"","RP","POI"};
  44
  45  // Standard magic:
  46  TString *baseDirPath = new TString(gSystem->pwd());
  47  TSystemDirectory *baseDir = new TSystemDirectory(".",baseDirPath->Data());
  48  TList *listOfFilesInBaseDir = baseDir->GetListOfFiles();
  49  TStopwatch timer;
  50  timer.Start();
  51  // listOfFilesInBaseDir->Print();
  52  Int_t nFiles = listOfFilesInBaseDir->GetEntries();
  53  Int_t fileCounter = 0;
  54  Double_t result[nMethods][3][nFilesMax] = {{{0.}}}; // [3 = "", "RP" or "POI"]
  55  Double_t error[nMethods][3][nFilesMax] = {{{0.}}}; // [3 = "", "RP" or "POI"]
  56  Double_t resultMinMax[nMethods][3][2] = {{{0.}}}; // [3 = "", "RP" or "POI"], [2 = min value, max value]
  57  for(Int_t m=0;m<nMethods;m++)
  58  {
  59   for(Int_t l=0;l<3;l++)
  60   {
  61    resultMinMax[m][l][0] = 44.;
  62    resultMinMax[m][l][1] = -44.;
  63   }
  64  }
  65  Double_t styleHistMinMax[3][2] = {{0.}}; // [3 = "", "RP" or "POI"], [2 = min value, max value]
  66  for(Int_t l=0;l<3;l++)
  67  {
  68   styleHistMinMax[l][0] = 44.;
  69   styleHistMinMax[l][1] = -44.;
  70  }
  71  cout<<endl;
  72  for(Int_t iFile=0;iFile<nFiles;iFile++)
  73  {
  74   TSystemFile *currentFile = (TSystemFile*)listOfFilesInBaseDir->At(iFile);
  75   // Consider only subdirectories:
  76   if(!currentFile ||
  77      !currentFile->IsDirectory() ||
  78      strcmp(currentFile->GetName(), ".") == 0 ||
  79     strcmp(currentFile->GetName(), "..") == 0) continue;
  80   // Accessing the output file "AnalysisResults.root" in current subdirectory:
  81   TString currentSubDirName = baseDirPath->Data();
  82   (currentSubDirName+="/")+=currentFile->GetName();
  83   currentSubDirName+="/";
  84   TString fileName = currentSubDirName;
  85   fileName+=outputFileName.Data();
  86   if(!(gSystem->AccessPathName(fileName.Data(),kFileExists)))
  87   {
  88    TFile *file = NULL;
  89    file = TFile::Open(fileName.Data(),"READ");
  90    for(Int_t m=0;m<nMethods;m++)
  91    {
  92     // Access from common output file the output file for particular method:
  93     TDirectoryFile *methodFile = NULL;
  94     methodFile = GetMethodFile(file,method[m],type);
  95     for(Int_t l=0;l<3;l++)
  96     {
  97      TH1D *histResult = NULL;
  98      histResult = GetHistogramWithResult(method[m],label[l],methodFile);
  99      if(histResult)
 100      {
 101       // Access the results:
 102       result[m][l][fileCounter] = histResult->GetBinContent(1);
 103       // Access the errors:
 104       error[m][l][fileCounter] = histResult->GetBinError(1);
 105       if(TMath::Abs(result[m][l][fileCounter])>pow(10.,-6.)) // take into account only if != 0 (to be improved - special care for < 0 is required)
 106       {
 107        // Establish min and max values for results:
 108        if(resultMinMax[m][l][0] > result[m][l][fileCounter]) resultMinMax[m][l][0] = result[m][l][fileCounter]; // min value
 109        if(resultMinMax[m][l][1] < result[m][l][fileCounter]) resultMinMax[m][l][1] = result[m][l][fileCounter]; // max value
 110        // Establish min and max values for style histograms:
 111        if(styleHistMinMax[l][0] > result[m][l][fileCounter]) styleHistMinMax[l][0] = result[m][l][fileCounter]; // min value
 112        if(styleHistMinMax[l][1] < result[m][l][fileCounter]) styleHistMinMax[l][1] = result[m][l][fileCounter]; // max value
 113       }
 114      }
 115     } // end of for(Int_t l=0;l<3;l++)
 116    } // end of for(Int_t m=0;m<nMethods;m++)
 117    if(file) file->Close();
 118    fileCounter++;
 119    //if(fileCounter%10==0)
 120    //{
 121     cout<<Form("Accessed %d files \"AnalysisResults.root\" so far....",fileCounter)<<"\r"<<flush;
 122    //}
 123   }
 124   if(fileCounter == nFilesMax) break;
 125  } // end of for(Int_t iFile=0;iFile<nFiles;iFile++)
 126
 127  cout<<Form("Accessed %d files \"AnalysisResults.root\" in total to estimate spread. ",fileCounter)<<endl;
 128  cout<<endl;
 129  if(fileCounter==0){exit(0);}
 130  const Int_t nFilesFinal = fileCounter;
 131
 132  // Make for each method graph holding results:
 133  TGraph *methodGraph[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"]
 134  Double_t x[nMethods][nFilesFinal] = {{0.}};
 135  for(Int_t m=0;m<nMethods;m++)
 136  {
 137   for(Int_t f=0;f<nFilesFinal;f++)
 138   {
 139    x[m][f]=m+0.5;
 140   }
 141  }
 142  for(Int_t m=0;m<nMethods;m++)
 143  {
 144   for(Int_t l=0;l<3;l++)
 145   {
 146    methodGraph[m][l] = new TGraph(nFilesFinal,x[m],result[m][l]);
 147    methodGraph[m][l]->SetMarkerStyle(methodMarkerStyle[m]);
 148    methodGraph[m][l]->SetMarkerColor(methodMarkerColor[m]);
 149   } // end of for(Int_t l=0;l<3;l++)
 150  } // for(Int_t m=0;m<nMethods;m++)
 151
 152  // Make for each method coloured mesh out of min and max values:
 153  Double_t meshWidth = 0.25;
 154  TGraph *methodMesh[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"]
 155  for(Int_t m=0;m<nMethods;m++)
 156  {
 157   for(Int_t l=0;l<3;l++)
 158   {
 159    if(resultMinMax[m][l][0]<44. && resultMinMax[m][l][1]>-44.)
 160    {
 161     methodMesh[m][l] = new TGraph(5);
 162     methodMesh[m][l]->SetPoint(0,(m+1-0.5)-meshWidth,resultMinMax[m][l][0]);
 163     methodMesh[m][l]->SetPoint(1,(m+1-0.5)+meshWidth,resultMinMax[m][l][0]);
 164     methodMesh[m][l]->SetPoint(2,(m+1-0.5)+meshWidth,resultMinMax[m][l][1]);
 165     methodMesh[m][l]->SetPoint(3,(m+1-0.5)-meshWidth,resultMinMax[m][l][1]);
 166     methodMesh[m][l]->SetPoint(4,(m+1-0.5)-meshWidth,resultMinMax[m][l][0]);
 167     methodMesh[m][l]->SetFillStyle(1001);
 168     methodMesh[m][l]->SetFillColor(methodMeshColor[m]);
 169    }
 170   }
 171  }
 172
 173  // Access for each method the results from the merged, large statistics file:
 174  Double_t resultMerged[nMethods][3] = {{0.}}; // [3 = "", "RP" or "POI"]
 175  TString mergedFileName = Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data());
 176  TFile *mergedFile = NULL;
 177  mergedFile = TFile::Open(mergedFileName.Data(),"READ");
 178  for(Int_t m=0;m<nMethods;m++)
 179  {
 180   TDirectoryFile *methodFile = NULL;
 181   if(!(gSystem->AccessPathName(fileName.Data(),kFileExists)))
 182   {
 183    if(mergedFile) methodFile = GetMethodFile(mergedFile,method[m],type);
 184   } else
 185     {
 186      cout<<"WARNING: Couldn't find the merged, large statistics file "<<endl;
 187      cout<<"         "<<fileName.Data()<<endl;
 188      cout<<"         in directory "<<gSystem->pwd()<<" !!!!"<<endl;
 189      cout<<"         Use macros mergeOuput.C and redoFinish.C to get it."<<endl;
 190      cout<<endl;
 191      break;
 192     }
 193   for(Int_t l=0;l<3;l++)
 194   {
 195    TH1D *histResult = NULL;
 196    if(methodFile)
 197    {
 198     histResult = GetHistogramWithResult(method[m],label[l],methodFile);
 199    }
 200    if(histResult)
 201    {
 202     // Access the results from the merged, large statistics file:
 203     resultMerged[m][l] = histResult->GetBinContent(1);
 204    } // end of for(Int_t l=0;l<3;l++)
 205   }
 206  } // end of for(Int_t m=0;m<nMethods;m++)
 207  if(mergedFile) mergedFile->Close();
 208
 209  // Make for each method graph holding results from the merged, large statistics file
 210  // and the errors from the randomly chosen small statistics file:
 211  TGraphErrors *methodMergedGraph[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"]
 212  Double_t xMerged[nMethods] = {0.};
 213  for(Int_t m=0;m<nMethods;m++)
 214  {
 215   xMerged[m]=m+0.5;
 216  }
 217  // Select randomly small statistics file:
 218  gRandom->SetSeed((UInt_t) (4400*timer.RealTime()/fileCounter));
 219  Int_t randomFile = (Int_t)gRandom->Uniform(0,fileCounter);
 220  for(Int_t m=0;m<nMethods;m++)
 221  {
 222   for(Int_t l=0;l<3;l++)
 223   {
 224    methodMergedGraph[m][l] = new TGraphErrors(1);
 225    methodMergedGraph[m][l]->SetPoint(0,xMerged[m],resultMerged[m][l]);
 226    if(showErrorOnMergedResult) methodMergedGraph[m][l]->SetPointError(0,0.,error[m][l][randomFile]);
 227    methodMergedGraph[m][l]->SetMarkerStyle(25);
 228    methodMergedGraph[m][l]->SetMarkerColor(kBlack);
 229   } // end of for(Int_t l=0;l<3;l++)
 230  } // for(Int_t m=0;m<nMethods;m++)
 231
 232  // Final drawing:
 233  gROOT->SetStyle("Plain"); // removing default gray color and setting white instead
 234  gStyle->SetOptStat(0); // removing statistics box from all histograms
 235  Bool_t showPlot[3] = {showPlotForReferenceFlow,showPlotForIntegratedFlowRP,showPlotForIntegratedFlowPOI};
 236  TString title[3] = {"Reference Flow","Integrated Flow (RP)","Integrated Flow (POI)"}
 237  TCanvas *canvas[3] = {NULL}; // [3 = "", "RP" or "POI"]
 238  for(Int_t l=0;l<3;l++)
 239  {
 240   if(showPlot[l])
 241   {
 242    canvas[l] = new TCanvas(Form("%s",title[l].Data()),Form("%s",title[l].Data()));
 243    TH1D *styleHist = StyleHist(title[l]);
 244    styleHist->SetMinimum(0.99*styleHistMinMax[l][0]);
 245    styleHist->SetMaximum(1.01*styleHistMinMax[l][1]);
 246    styleHist->Draw();
 247    for(Int_t m=0;m<nMethods;m++)
 248    {
 249     if(methodMesh[m][l]) methodMesh[m][l]->Draw("lfsame");
 250     if(methodGraph[m][l]) methodGraph[m][l]->Draw("psame");
 251     if(TMath::Abs(*(methodMergedGraph[m][l]->GetY()))>pow(10.,-6.)) // draw only if not == 0.
 252     {
 253      methodMergedGraph[m][l]->Draw("psame");
 254     }
 255    } // end of for(Int_t m=0;m<nMethods;m++)
 256   } // end of if(showPlot[l])
 257  } // end of for(Int_t l=0;l<3;l++)
 258
 259  timer.Stop();
 260  timer.Print();
 261  cout<<endl;
 262
 263 } // end of void showSpread(TString type="", Int_t mode=mLocal)
 264
 265 // =============================================================================================
 266
 267 TH1D *StyleHist(TString histTitle)
 268 {
 269  // Make the style histogram.
 270  Int_t n = 2; // harmonic (to be improved - access this from common control hist)
 271  TH1D *styleHist = new TH1D(Form("%s",histTitle.Data()),Form("%s",histTitle.Data()),nMethods,0,nMethods);
 272  styleHist->GetYaxis()->SetTickLength(0.01);
 273  for(Int_t m=0;m<nMethods;m++)
 274  {
 275   styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,method[m].Data()));
 276   if(method[m]=="LYZ1SUM" || method[m]=="LYZ2SUM")
 277   {
 278    styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,"LYZ,sum"));
 279   }
 280   else if(method[m]=="LYZ1PROD" || method[m]=="LYZ2PROD")
 281   {
 282    styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,"LYZ,prod"));
 283   }
 284  } // end of for(Int_t m=0;m<nMethods;m++)
 285
 286  return styleHist;
 287 }
 288
 289 // =============================================================================================
 290
 291 TDirectoryFile* GetMethodFile(TFile *commonFile, TString method, TString type)
 292 {
 293  // Form a file name for each method:
 294  TString methodFileName = "output";
 295  if(method.Contains("GFC"))
 296  {
 297   methodFileName+="GFC";
 298  } else if(method.Contains("QC"))
 299    {
 300     methodFileName+="QC";
 301    } else
 302      {
 303       methodFileName+=method.Data();
 304      }
 305  methodFileName+="analysis";
 306  methodFileName+=type.Data();
 307  TDirectoryFile *methodFile = NULL;
 308  if(commonFile)
 309  {
 310   methodFile = (TDirectoryFile*)commonFile->FindObjectAny(methodFileName.Data());
 311  }
 312
 313  return methodFile;
 314
 315 } // end of TDirectoryFile* AccessMethodFile(TString commonFile, TString method, TString type)
 316
 317 // =============================================================================================
 318
 319 TH1D* GetHistogramWithResult(TString method, TString label, TDirectoryFile *methodFile)
 320 {
 321  // Access first the common list holding all output histograms:
 322  TList *methodList = NULL;
 323  if(method.Contains("GFC"))
 324  {
 325   methodFile->GetObject("cobjGFC",methodList);
 326  } else if(method.Contains("QC"))
 327    {
 328     methodFile->GetObject("cobjQC",methodList);
 329    } else
 330      {
 331       methodFile->GetObject(Form("cobj%s",method.Data()),methodList);
 332      }
 333  // Access from the common list the needed histogram:
 334  if(methodList)
 335  {
 336   AliFlowCommonHistResults *commonHistRes = NULL;
 337   if(!(method.Contains("GFC") || method.Contains("QC")))
 338   {
 339    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%s",method.Data()));
 340   }
 341   else if(method=="2,GFC")
 342   {
 343    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults2ndOrderGFC");
 344   }
 345   else if(method=="4,GFC")
 346   {
 347    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults4thOrderGFC");
 348   }
 349   else if(method=="6,GFC")
 350   {
 351    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults6thOrderGFC");
 352   }
 353   else if(method=="8,GFC")
 354   {
 355    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults8thOrderGFC");
 356   }
 357   else if(method=="2,QC")
 358   {
 359    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults2ndOrderQC");
 360   }
 361   else if(method=="4,QC")
 362   {
 363    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults4thOrderQC");
 364   }
 365   else if(method=="6,QC")
 366   {
 367    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults6thOrderQC");
 368   }
 369   else if(method=="8,QC")
 370   {
 371    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults8thOrderQC");
 372   }
 373  } // end of if(methodList)
 374
 375  // Access histogram with results for reference flow or integrated flow of RPs or POIs:
 376  TH1D *hist = NULL;
 377  if(label=="")
 378  {
 379   hist = commonHistRes->GetHistIntFlow();
 380  } else if(label=="RP")
 381    {
 382     hist = commonHistRes->GetHistIntFlowRP();
 383    } else if(label=="POI")
 384      {
 385       hist = commonHistRes->GetHistIntFlowPOI();
 386      }
 387
 388  if(hist) return hist;
 389
 390 } // end of TH1D *GetHistogramWithResult(TString method, TString rf_rp_poi, TString file);
 391
 392 // =============================================================================================
 393
 394 void CrossCheckUserSettings()
 395 {
 396  // Check in this method if the user settings make sense.
 397
 398  if(nFilesMax<=0)
 399  {
 400   cout<<endl;
 401   cout<<"WARNING: nFilesMax must be a positive integer (not too large, though) !!!!"<<endl;
 402   cout<<endl;
 403   exit(0);
 404  }
 405  if(nFilesMax>44)
 406  {
 407   cout<<endl;
 408   cout<<"WARNING: You may want to set nFilesMax to the smaller value."<<endl;
 409   cout<<"         Otherwise you might wait forever to see the plots."<<endl;
 410  }
 411
 412 } // end of void CrossCheckUserSettings()
 413
 414 // =============================================================================================
 415
 416 void LoadLibrariesSS(const libModes mode) {
 417
 418   //--------------------------------------
 419   // Load the needed libraries most of them already loaded by aliroot
 420   //--------------------------------------
 421   //gSystem->Load("libTree");
 422   gSystem->Load("libGeom");
 423   gSystem->Load("libVMC");
 424   gSystem->Load("libXMLIO");
 425   gSystem->Load("libPhysics");
 426
 427   //----------------------------------------------------------
 428   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 429   //----------------------------------------------------------
 430   if (mode==mLocal) {
 431     //--------------------------------------------------------
 432     // If you want to use already compiled libraries
 433     // in the aliroot distribution
 434     //--------------------------------------------------------
 435
 436   //==================================================================================
 437   //load needed libraries:
 438   gSystem->AddIncludePath("-I$ROOTSYS/include");
 439   //gSystem->Load("libTree");
 440
 441   // for AliRoot
 442   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 443   gSystem->Load("libANALYSIS");
 444   gSystem->Load("libPWGflowBase");
 445   //cerr<<"libPWGflowBase loaded ..."<<endl;
 446
 447   }
 448
 449   else if (mode==mLocalSource) {
 450
 451     // In root inline compile
 452
 453     // Constants
 454     gROOT->LoadMacro("Base/AliFlowCommonConstants.cxx+");
 455     gROOT->LoadMacro("Base/AliFlowLYZConstants.cxx+");
 456
 457     // Flow event
 458     gROOT->LoadMacro("Base/AliFlowVector.cxx+");
 459     gROOT->LoadMacro("Base/AliFlowTrackSimple.cxx+");
 460     gROOT->LoadMacro("Base/AliFlowTrackSimpleCuts.cxx+");
 461     gROOT->LoadMacro("Base/AliFlowEventSimple.cxx+");
 462
 463     // Output histograms
 464     gROOT->LoadMacro("Base/AliFlowCommonHist.cxx+");
 465     gROOT->LoadMacro("Base/AliFlowCommonHistResults.cxx+");
 466     gROOT->LoadMacro("Base/AliFlowLYZHist1.cxx+");
 467     gROOT->LoadMacro("Base/AliFlowLYZHist2.cxx+");
 468
 469     // Functions needed for various methods
 470     gROOT->LoadMacro("Base/AliCumulantsFunctions.cxx+");
 471     gROOT->LoadMacro("Base/AliFlowLYZEventPlane.cxx+");
 472
 473     // Flow Analysis code for various methods
 474     gROOT->LoadMacro("Base/AliFlowAnalysisWithMCEventPlane.cxx+");
 475     gROOT->LoadMacro("Base/AliFlowAnalysisWithScalarProduct.cxx+");
 476     gROOT->LoadMacro("Base/AliFlowAnalysisWithLYZEventPlane.cxx+");
 477     gROOT->LoadMacro("Base/AliFlowAnalysisWithLeeYangZeros.cxx+");
 478     gROOT->LoadMacro("Base/AliFlowAnalysisWithCumulants.cxx+");
 479     gROOT->LoadMacro("Base/AliFlowAnalysisWithQCumulants.cxx+");
 480     gROOT->LoadMacro("Base/AliFlowAnalysisWithFittingQDistribution.cxx+");
 481     gROOT->LoadMacro("Base/AliFlowAnalysisWithMixedHarmonics.cxx+");
 482     gROOT->LoadMacro("Base/AliFlowAnalysisWithNestedLoops.cxx+");
 483
 484     cout << "finished loading macros!" << endl;
 485
 486   } // end of else if (mode==mLocalSource)
 487
 488 } // end of void LoadLibrariesSS(const libModes mode)