ROOT logo
//==========================================================================//
// Macro to refit the q-distribution in FQD method. By using this macro the // 
// fitting parameters for q-distribution can be tuned and correspondingly   //
// the estimate for reference flow harmonic obtained from FQD method can be //
// significantly improved.                                                  // 
//==========================================================================//

// Fitting parameters for q-distribution:
Double_t vStart = 0.05; // starting value for v fit
Double_t vMin = 0.0; // lower bound for v fit 
Double_t vMax = 0.25; // upper bound for v fit  
Double_t sigma2Start = 0.75; // starting value for sigma^2 fit 
Double_t sigma2Min = 0.5; // lower bound for sigma^2 fit (according to theorists must be >= 0.5)   
Double_t sigma2Max = 2.5; // upper bound for sigma^2 fit
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
Bool_t finalResultIsFromSigma2Fitted = kTRUE; // final saved result is obtained with sigma^2 fitted (kTRUE) or sigma^2 fixed (kFALSE)
Bool_t printOnTheScreen = kTRUE; // print or not the final results on the screen
Bool_t plotResults = kTRUE; // plot on the screen q-distribution and resulting fitting functions (for non-fixed sigma^2 and fixed sigma^2)

// Name of the common output file:
TString outputFileName = "AnalysisResults.root";
//TString outputFileName = "outputCentrality0.root";

enum libModes {mLocal,mLocalSource};

TFile *commonOutputFile = NULL; // common output file "AnalysisResults.root"
TDirectoryFile *dirFileFQD = NULL; // FQD's TDirectoryFile in "AnalysisResults.root"
TList *outputListFQD = NULL; // output list holding all FQD objects

void fqd(TString analysisType="ESD", Int_t analysisMode=mLocal)
{
 // 1. analysisType: type of analysis can be ESD, AOD, MC, MK, ....
 //                  (if type="" output files are from simulation 'on the fly')
 // 2. analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
 //                  if analysisMode = mLocalSource -> analyze data on your computer using root + source files 
  
 // Cross-check if the user's settings make sense:
 CrossCheckSettings();
 // Load needed libraries:                       
 LoadLibrariesFQD(analysisMode); 
 // Get output list which holds all FQD objects:
 GetOutputList(analysisType); 
 // Redo fit of q-distribution:
 RedoFit();
 // Plot q-distribution and resulting fitting functions:
 if(plotResults) Plot(); 
 // Save the new results in the common output file: 
 Save();
  
} // end of void fqd(TString type="ESD", Int_t mode=mLocal)   

// =========================================================================================== 

void RedoFit()
{
 // Redo the fit of q-distribution.

 AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 fqd->GetOutputHistograms(outputListFQD);
 // Set new fitting parameters:
 fqd->SetvStart(vStart);  
 fqd->SetvMin(vMin);
 fqd->SetvMax(vMax);
 fqd->SetSigma2Start(sigma2Start);  
 fqd->SetSigma2Min(sigma2Min);
 fqd->SetSigma2Max(sigma2Max);
 fqd->SetTreshold(treshold);
 fqd->SetFinalResultIsFromSigma2Fitted(finalResultIsFromSigma2Fitted);
 fqd->SetPrintOnTheScreen(printOnTheScreen); 
 // Save new fitting parameters:  
 fqd->StoreFittingParameters();
 // Redo fit:
 fqd->Finish(kTRUE);
 
} // end of void RedoFit(TList *outputList)
 
// =========================================================================================== 
 
void GetOutputList(TString analysisType)
{
 // Get output list which holds all FQD objects.
  
 // Access common output file:
 commonOutputFile = AccessOutputFile(outputFileName);
 
 // Access from common output file the TDirectoryFile for FQD method
 // and from it the output list holding all objects:
 GetListWithHistograms(commonOutputFile,analysisType,"FQD");
   
} // end of TList* GetOutputList(TString analysisType)
 
// ===========================================================================================

void Plot()
{
 // Plot q-distribution and resulting fitting functions.
 
 gROOT->SetStyle("Plain"); // default color is white instead of gray
 gStyle->SetOptStat(0); // remove stat. box from all histos
 legend = new TLegend(0.6,0.55,0.85,0.7); 
 legend->SetFillStyle(0);
 // q-distribution:
 TH1D *qDistribution = dynamic_cast<TH1D*>(outputListFQD->FindObject("fqDistribution"));
 Cosmetics(qDistribution);
 legend->AddEntry(qDistribution,"q-distribution","f");
 qDistribution->Draw();
 // resulting fitting functions:
 TF1 *fittingFun[2] = {NULL};
 TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
 Double_t lineColors[2] = {kBlue,kRed};
 for(Int_t s2F=0;s2F<2;s2F++)
 {
  fittingFun[s2F] = dynamic_cast<TF1*>(outputListFQD->FindObject(Form("fFittingFunction, %s",sigmaFlag[s2F].Data())));
  fittingFun[s2F]->SetLineColor(lineColors[s2F]);
  fittingFun[s2F]->Draw("SAME");
  legend->AddEntry(fittingFun[s2F]->GetName(),sigmaFlag[s2F].Data(),"l");
 }
 legend->Draw("SAME"); 
  
} // end of Plot()

// ===========================================================================================

void Save()
{
 // Save the new results of fitting for FQD method in the common output file.
 
 dirFileFQD->Add(outputListFQD,kTRUE);
 dirFileFQD->Write(dirFileFQD->GetName(),TObject::kSingleKey+TObject::kOverwrite);
 //delete commonOutputFile;

} // end of void Save()

// ===========================================================================================

void Cosmetics(TH1D *hist)
{
 // Set cosmetics for the q-distribution.
 
 Int_t firstNonEmptyBin = hist->FindFirstBinAbove(0);
 Double_t lowRange = hist->GetBinLowEdge(firstNonEmptyBin);
 Int_t lastNonEmptyBin = hist->FindLastBinAbove(0);
 Double_t upperRange = hist->GetBinLowEdge(lastNonEmptyBin+10);
 hist->GetXaxis()->SetRangeUser(lowRange,upperRange);  
 hist->SetFillColor(16);  
 hist->SetTitle("Fitted q-distribution");
 
} // end of void Cosmetics(TH1D *hist)

// ===========================================================================================

TFile* AccessOutputFile(TString outputFileName)
{
 // Access the common output file.
    
 TFile *outputFile = NULL;
 if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
 {
  outputFile = TFile::Open(outputFileName.Data(),"UPDATE");
 } else
   {
    cout<<endl;
    cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl;
    cout<<"         directory "<<gSystem->pwd()<<" !!!!"<<endl;
    cout<<endl;
    exit(0);
   }
  
 return outputFile;
 
} // end of TFile* AccessOutputFile(TString outputFileName)
  
// ===========================================================================================
 
void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
{
 // Access from common output file the TDirectoryFile for FQD method and from it
 // the output list holding all objects:

 TString fileName = ""; 
 TString listName = ""; 
 // Form a file name: 
 fileName+="output";
 fileName+=methodName.Data();
 fileName+="analysis";
 fileName+=analysisType.Data();
 // Access this file:
 dirFileFQD = (TDirectoryFile*)outputFile->FindObjectAny(fileName.Data());
 // Form a list name for each method:
 if(dirFileFQD)
 {
  TList* listTemp = dirFileFQD->GetListOfKeys();
  if(listTemp && listTemp->GetEntries() == 1)
  {
   listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast)
   dirFileFQD->GetObject(listName.Data(),outputListFQD);
  } else
    {
     cout<<" WARNING: Accessing TList from TDirectoryFile failed miserably for FQD method !!!!"<<endl;
     cout<<"          Did you actually use FQD in the analysis?"<<endl;
     cout<<endl;
     exit(0);
    }
 } else 
   {
    cout<<" WARNING: Couldn't find a TDirectoryFile "<<fileName.Data()<<" in file "<<outputFileName.Data()<<" !!!!"<<endl;
    cout<<"          Perhaps wrong analysis type? Can be \"ESD\",\"AOD\",\"\", etc."<<endl;
    cout<<"          Or you didn't use FQD in the analysis which have produced "<<outputFileName.Data()<<"?"<<endl;
   } 
          
} // end of void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
 
//===========================================================================================

void CrossCheckSettings() 
{
 // Check in this method if the settings make sense.
 
 if(sigma2Min<0.5)
 {
  cout<<endl; 
  cout<<"WARNING: According to theorists sigma2Min >= 0.5 !!!!"<<endl;
  cout<<endl; 
  exit(0);
 }
 
} // end of void CrossCheckSettings()

//===========================================================================================

void LoadLibrariesFQD(const libModes mode) {
  
  //--------------------------------------
  // Load the needed libraries most of them already loaded by aliroot
  //--------------------------------------
  gSystem->Load("libTree.so");
  gSystem->Load("libGeom.so");
  gSystem->Load("libVMC.so");
  gSystem->Load("libXMLIO.so");
  gSystem->Load("libPhysics.so");
  
  //----------------------------------------------------------
  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
  //----------------------------------------------------------
  if (mode==mLocal) {
    //--------------------------------------------------------
    // If you want to use already compiled libraries 
    // in the aliroot distribution
    //--------------------------------------------------------

  //==================================================================================  
  //load needed libraries:
  gSystem->AddIncludePath("-I$ROOTSYS/include");
  gSystem->Load("libTree.so");

  // for AliRoot
  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libPWGflowBase.so");
  cerr<<"libPWGflowBase.so loaded ..."<<endl;
  
  }
  
  else if (mode==mLocalSource) {
 
    // In root inline compile
  
    // Constants  
    gROOT->LoadMacro("BaseAliFlowCommonConstants.cxx+");
    gROOT->LoadMacro("BaseAliFlowLYZConstants.cxx+");
    
    // Flow event
    gROOT->LoadMacro("BaseAliFlowVector.cxx+"); 
    gROOT->LoadMacro("BaseAliFlowTrackSimple.cxx+");    
    gROOT->LoadMacro("BaseAliFlowTrackSimpleCuts.cxx+");    
    gROOT->LoadMacro("BaseAliFlowEventSimple.cxx+");
    
    // Output histosgrams
    gROOT->LoadMacro("BaseAliFlowCommonHist.cxx+");
    gROOT->LoadMacro("BaseAliFlowCommonHistResults.cxx+");
    gROOT->LoadMacro("BaseAliFlowLYZHist1.cxx+");
    gROOT->LoadMacro("BaseAliFlowLYZHist2.cxx+");
       
    cout << "finished loading macros!" << endl;  
    
  } // end of else if (mode==mLocalSource) 
  
} // end of void LoadLibrariesFQD(const libModes mode)
 fqd.C:1
 fqd.C:2
 fqd.C:3
 fqd.C:4
 fqd.C:5
 fqd.C:6
 fqd.C:7
 fqd.C:8
 fqd.C:9
 fqd.C:10
 fqd.C:11
 fqd.C:12
 fqd.C:13
 fqd.C:14
 fqd.C:15
 fqd.C:16
 fqd.C:17
 fqd.C:18
 fqd.C:19
 fqd.C:20
 fqd.C:21
 fqd.C:22
 fqd.C:23
 fqd.C:24
 fqd.C:25
 fqd.C:26
 fqd.C:27
 fqd.C:28
 fqd.C:29
 fqd.C:30
 fqd.C:31
 fqd.C:32
 fqd.C:33
 fqd.C:34
 fqd.C:35
 fqd.C:36
 fqd.C:37
 fqd.C:38
 fqd.C:39
 fqd.C:40
 fqd.C:41
 fqd.C:42
 fqd.C:43
 fqd.C:44
 fqd.C:45
 fqd.C:46
 fqd.C:47
 fqd.C:48
 fqd.C:49
 fqd.C:50
 fqd.C:51
 fqd.C:52
 fqd.C:53
 fqd.C:54
 fqd.C:55
 fqd.C:56
 fqd.C:57
 fqd.C:58
 fqd.C:59
 fqd.C:60
 fqd.C:61
 fqd.C:62
 fqd.C:63
 fqd.C:64
 fqd.C:65
 fqd.C:66
 fqd.C:67
 fqd.C:68
 fqd.C:69
 fqd.C:70
 fqd.C:71
 fqd.C:72
 fqd.C:73
 fqd.C:74
 fqd.C:75
 fqd.C:76
 fqd.C:77
 fqd.C:78
 fqd.C:79
 fqd.C:80
 fqd.C:81
 fqd.C:82
 fqd.C:83
 fqd.C:84
 fqd.C:85
 fqd.C:86
 fqd.C:87
 fqd.C:88
 fqd.C:89
 fqd.C:90
 fqd.C:91
 fqd.C:92
 fqd.C:93
 fqd.C:94
 fqd.C:95
 fqd.C:96
 fqd.C:97
 fqd.C:98
 fqd.C:99
 fqd.C:100
 fqd.C:101
 fqd.C:102
 fqd.C:103
 fqd.C:104
 fqd.C:105
 fqd.C:106
 fqd.C:107
 fqd.C:108
 fqd.C:109
 fqd.C:110
 fqd.C:111
 fqd.C:112
 fqd.C:113
 fqd.C:114
 fqd.C:115
 fqd.C:116
 fqd.C:117
 fqd.C:118
 fqd.C:119
 fqd.C:120
 fqd.C:121
 fqd.C:122
 fqd.C:123
 fqd.C:124
 fqd.C:125
 fqd.C:126
 fqd.C:127
 fqd.C:128
 fqd.C:129
 fqd.C:130
 fqd.C:131
 fqd.C:132
 fqd.C:133
 fqd.C:134
 fqd.C:135
 fqd.C:136
 fqd.C:137
 fqd.C:138
 fqd.C:139
 fqd.C:140
 fqd.C:141
 fqd.C:142
 fqd.C:143
 fqd.C:144
 fqd.C:145
 fqd.C:146
 fqd.C:147
 fqd.C:148
 fqd.C:149
 fqd.C:150
 fqd.C:151
 fqd.C:152
 fqd.C:153
 fqd.C:154
 fqd.C:155
 fqd.C:156
 fqd.C:157
 fqd.C:158
 fqd.C:159
 fqd.C:160
 fqd.C:161
 fqd.C:162
 fqd.C:163
 fqd.C:164
 fqd.C:165
 fqd.C:166
 fqd.C:167
 fqd.C:168
 fqd.C:169
 fqd.C:170
 fqd.C:171
 fqd.C:172
 fqd.C:173
 fqd.C:174
 fqd.C:175
 fqd.C:176
 fqd.C:177
 fqd.C:178
 fqd.C:179
 fqd.C:180
 fqd.C:181
 fqd.C:182
 fqd.C:183
 fqd.C:184
 fqd.C:185
 fqd.C:186
 fqd.C:187
 fqd.C:188
 fqd.C:189
 fqd.C:190
 fqd.C:191
 fqd.C:192
 fqd.C:193
 fqd.C:194
 fqd.C:195
 fqd.C:196
 fqd.C:197
 fqd.C:198
 fqd.C:199
 fqd.C:200
 fqd.C:201
 fqd.C:202
 fqd.C:203
 fqd.C:204
 fqd.C:205
 fqd.C:206
 fqd.C:207
 fqd.C:208
 fqd.C:209
 fqd.C:210
 fqd.C:211
 fqd.C:212
 fqd.C:213
 fqd.C:214
 fqd.C:215
 fqd.C:216
 fqd.C:217
 fqd.C:218
 fqd.C:219
 fqd.C:220
 fqd.C:221
 fqd.C:222
 fqd.C:223
 fqd.C:224
 fqd.C:225
 fqd.C:226
 fqd.C:227
 fqd.C:228
 fqd.C:229
 fqd.C:230
 fqd.C:231
 fqd.C:232
 fqd.C:233
 fqd.C:234
 fqd.C:235
 fqd.C:236
 fqd.C:237
 fqd.C:238
 fqd.C:239
 fqd.C:240
 fqd.C:241
 fqd.C:242
 fqd.C:243
 fqd.C:244
 fqd.C:245
 fqd.C:246
 fqd.C:247
 fqd.C:248
 fqd.C:249
 fqd.C:250
 fqd.C:251
 fqd.C:252
 fqd.C:253
 fqd.C:254
 fqd.C:255
 fqd.C:256
 fqd.C:257
 fqd.C:258
 fqd.C:259
 fqd.C:260
 fqd.C:261
 fqd.C:262
 fqd.C:263
 fqd.C:264
 fqd.C:265
 fqd.C:266
 fqd.C:267
 fqd.C:268
 fqd.C:269
 fqd.C:270
 fqd.C:271
 fqd.C:272
 fqd.C:273
 fqd.C:274
 fqd.C:275
 fqd.C:276
 fqd.C:277
 fqd.C:278
 fqd.C:279
 fqd.C:280
 fqd.C:281
 fqd.C:282
 fqd.C:283
 fqd.C:284
 fqd.C:285
 fqd.C:286
 fqd.C:287
 fqd.C:288
 fqd.C:289