PWGCF/FLOW/macros/parity/plotMH.C

   1 //=====================================================//
   2 //Macro that reads the output of the flow analysis and
   3 //retrieves the QC and the MH containers.
   4 //The macro produces an output called outputMH.root
   5 //that contains the four histograms related to the
   6 //3 particle correlator and the two integrated v2{2,QC}
   7 //and v2{4,QC}
   8 //Author: Panos.Christakoglou@cern.ch
   9
  10 void plotMH(const char* filename = "AnalysisResults.root",
  11             const char* analysisType = "TPC only",
  12             Int_t centrality = 0) {
  13   gStyle->SetPalette(1,0);
  14
  15   //----------------------------------------------------------
  16   // >>>>>>>>>>> Load libraries <<<<<<<<<<<<<<
  17   //----------------------------------------------------------
  18   gSystem->AddIncludePath("-I$ROOTSYS/include");
  19   gSystem->Load("libGeom");
  20   gSystem->Load("libVMC");
  21   gSystem->Load("libXMLIO");
  22   gSystem->Load("libPhysics");
  23
  24   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  25   gSystem->Load("libANALYSIS");
  26   gSystem->Load("libPWGflowBase");
  27
  28   //----------------------------------------------------------
  29   // >>>>>>>>>>> Open file - Get objects <<<<<<<<<<<<<<
  30   //----------------------------------------------------------
  31   TFile *fInput = TFile::Open(filename);
  32   if(!fInput) {
  33     Printf("File not found!!!");
  34     break;
  35   }
  36   //f->ls();
  37
  38   //Get the TList of the MH
  39   TList *listMH = GetMHResults(fInput,analysisType,centrality);
  40
  41   //Get the TList of the QC
  42   TList *listQC = GetQCResults(fInput,analysisType,centrality);
  43
  44   TFile *fOutput = TFile::Open("outputMH.root","recreate");
  45   listMH->Write();
  46   listQC->Write();
  47   fOutput->Close();
  48 }
  49
  50 //____________________________________________________________//
  51 TList *GetQCResults(TFile *fInput,
  52                     const char* analysisType = 0x0,
  53                     Int_t centrality = -1) {
  54   //Function that reads the TDirectoryFile of the MH
  55   //and returns a TList with the relevant plots.
  56   TList *listOutput = new TList();
  57   listOutput->SetName("listQC");
  58
  59   //______________________________________________________________//
  60   //Global variables
  61   const Int_t nQCMethods = 4;
  62   AliFlowCommonHistResults *commonHistRes[nQCMethods];
  63   TString methodIntegratedFlowQC[nQCMethods] = {"2ndOrderQC",
  64                                                 "4thOrderQC",
  65                                                 "6thOrderQC",
  66                                                 "8thOrderQC"};
  67   TString gAliFlowCommonHistName[nQCMethods];
  68
  69   //______________________________________________________________//
  70   //Get the TDirectoryFile
  71   TString directoryNameQC = 0;
  72   directoryNameQC = "outputQCanalysis";
  73   if(analysisType) directoryNameQC += analysisType;
  74   TDirectoryFile *outputQCanalysis = dynamic_cast<TDirectoryFile *>(fInput->Get(directoryNameQC.Data()));
  75   if(!outputQCanalysis) {
  76     Printf("QC directory not found!!!");
  77     break;
  78   }
  79   //outputQCanalysis->ls();
  80
  81   //______________________________________________________________//
  82   //Get the TList
  83   TString listNameQC = 0;
  84   if(centrality > -1) {
  85     listNameQC = "cobjQC_outputCentrality";
  86     listNameQC += centrality; listNameQC += ".root";
  87   }
  88   else listNameQC = "cobjQC";
  89   TList *cobjQC = dynamic_cast<TList *>(outputQCanalysis->Get(listNameQC.Data()));
  90   if(!cobjQC) {
  91     Printf("QC object list not found!!!");
  92     break;
  93   }
  94
  95   //Common hist for QC
  96   Double_t nAnalyzedEvents = 0.;
  97   AliFlowCommonHist *commonHistQC = dynamic_cast<AliFlowCommonHist *>(cobjQC->FindObject("AliFlowCommonHistQC"));
  98   if(commonHistQC) {
  99     TList *clistQCCommonHist = dynamic_cast<TList *>(commonHistQC->GetHistList());
 100     if(clistQCCommonHist) {
 101       TH1F *gHistMultRPQC = dynamic_cast<TH1F *>(clistQCCommonHist->FindObject("Control_Flow_MultRP_AliFlowCommonHistQC"));
 102       if(gHistMultRPQC)
 103         nAnalyzedEvents = gHistMultRPQC->GetEntries();
 104     }
 105   }
 106
 107   //2nd order QC
 108   gAliFlowCommonHistName[0] = "AliFlowCommonHistResults";
 109   gAliFlowCommonHistName[0] += methodIntegratedFlowQC[0].Data();
 110   commonHistRes[0] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[0].Data()));
 111   //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[0].Data());
 112   TH1D *histQC_2 = commonHistRes[0]->GetHistIntFlowRP();
 113   histQC_2->SetName("fHistIntegratedFlowRPQC_2");
 114
 115   //4th order QC
 116   gAliFlowCommonHistName[1] = "AliFlowCommonHistResults";
 117   gAliFlowCommonHistName[1] += methodIntegratedFlowQC[1].Data();
 118   commonHistRes[1] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[1].Data()));
 119   //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[1].Data());
 120   TH1D *histQC_4 = commonHistRes[1]->GetHistIntFlowRP();
 121   histQC_4->SetName("fHistIntegratedFlowRPQC_4");
 122
 123   //6th order QC
 124   gAliFlowCommonHistName[2] = "AliFlowCommonHistResults";
 125   gAliFlowCommonHistName[2] += methodIntegratedFlowQC[2].Data();
 126   commonHistRes[2] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[2].Data()));
 127   //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[2].Data());
 128   TH1D *histQC_6 = commonHistRes[2]->GetHistIntFlowRP();
 129   histQC_6->SetName("fHistIntegratedFlowRPQC_6");
 130
 131   //8th order QC
 132   gAliFlowCommonHistName[3] = "AliFlowCommonHistResults";
 133   gAliFlowCommonHistName[3] += methodIntegratedFlowQC[3].Data();
 134   commonHistRes[3] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[3].Data()));
 135   //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[3].Data());
 136   TH1D *histQC_8 = commonHistRes[3]->GetHistIntFlowRP();
 137   histQC_8->SetName("fHistIntegratedFlowRPQC_8");
 138
 139   TH1D *gHistEvents = new TH1D("gHistEvents",";;N_{analyzed events}",
 140                                1,0.5,1.5);
 141   gHistEvents->SetBinContent(1,nAnalyzedEvents);
 142   //______________________________________________________________//
 143   //Print the results
 144   Printf("\n============================================================");
 145   Printf("Analyzed events: %d",(Int_t)nAnalyzedEvents);
 146   Printf("v2(2,QC): %lf +- %lf",histQC_2->GetBinContent(1),
 147          histQC_2->GetBinError(1));
 148   Printf("v2(4,QC): %lf +- %lf",histQC_4->GetBinContent(1),
 149          histQC_4->GetBinError(1));
 150   Printf("v2(6,QC): %lf +- %lf",histQC_6->GetBinContent(1),
 151          histQC_6->GetBinError(1));
 152   Printf("v2(8,QC): %lf +- %lf",histQC_8->GetBinContent(1),
 153          histQC_8->GetBinError(1));
 154   Printf("============================================================");
 155
 156   listOutput->Add(histQC_2);
 157   listOutput->Add(histQC_4);
 158   listOutput->Add(histQC_6);
 159   listOutput->Add(histQC_8);
 160   listOutput->Add(gHistEvents);
 161
 162   return listOutput;
 163 }
 164
 165 //____________________________________________________________//
 166 TList *GetMHResults(TFile *fInput,
 167                     const char* analysisType = 0x0,
 168                     Int_t centrality = -1) {
 169   //Function that reads the TDirectoryFile of the MH
 170   //and returns a TList with the relevant plots.
 171   TList *listOutput = new TList();
 172   listOutput->SetName("listMH");
 173
 174   //______________________________________________________________//
 175   //Get the TDirectoryFile
 176   TString directoryNameMH = 0;
 177   directoryNameMH = "outputMHanalysis";
 178   if(analysisType) directoryNameMH += analysisType;
 179   TDirectoryFile *outputMHanalysis = dynamic_cast<TDirectoryFile *>(fInput->Get(directoryNameMH.Data()));
 180   if(!outputMHanalysis) {
 181     Printf("MH directory not found!!!");
 182     break;
 183   }
 184   //outputMHanalysis->ls();
 185
 186   //______________________________________________________________//
 187   //Get the TList
 188   TString listNameMH = 0;
 189   if(centrality > -1) {
 190     listNameMH = "cobjMH_outputCentrality";
 191     listNameMH += centrality; listNameMH += ".root";
 192   }
 193   else listNameMH = "cobjMH";
 194   TList *cobjMH = dynamic_cast<TList *>(outputMHanalysis->Get(listNameMH.Data()));
 195   if(!cobjMH) {
 196     Printf("MH object list not found!!!");
 197     break;
 198   }
 199   //cobjMH->ls();
 200
 201   //______________________________________________________________//
 202   //Get the daughter TLists
 203   TList *listWeights = dynamic_cast<TList *>(cobjMH->At(0));
 204   //listWeights->ls();
 205   TList *listProfiles = dynamic_cast<TList *>(cobjMH->At(1));
 206   //listProfiles->ls();
 207   TList *listResults = dynamic_cast<TList *>(cobjMH->At(2));
 208   //listResults->ls();
 209
 210   if((!listWeights)||(!listProfiles)||(!listResults)) {
 211     Printf("MH output lists not found!!!");
 212     break;
 213   }
 214
 215   //______________________________________________________________//
 216   TProfile *fAnalysisSettings = dynamic_cast<TProfile *>(cobjMH->At(3));
 217
 218   //______________________________________________________________//
 219   //Get the objects from the Results list
 220   TH1D *f3pCorrelatorHist = dynamic_cast<TH1D *>(listResults->At(0));
 221   TH1D *fDetectorBiasHist = dynamic_cast<TH1D *>(listResults->At(1));
 222   TH1D *f3pCorrelatorVsMHist = dynamic_cast<TH1D *>(listResults->At(2));
 223   TH1D *fDetectorBiasVsMHist = dynamic_cast<TH1D *>(listResults->At(3));
 224
 225   //______________________________________________________________//
 226   //Get the objects from the Profile list
 227   TProfile *f3pCorrelatorPro = dynamic_cast<TProfile *>(listProfiles->At(0));
 228   TProfile *fNonIsotropicTermsPro = dynamic_cast<TProfile *>(listProfiles->At(1));
 229   TProfile *f3pCorrelatorVsMPro = dynamic_cast<TProfile *>(listProfiles->At(2));
 230   TProfile2D *fNonIsotropicTermsVsMPro = dynamic_cast<TProfile2D *>(listProfiles->At(3));
 231   TProfile *f3pCorrelatorVsPtSumPro = dynamic_cast<TProfile *>(listProfiles->At(4));
 232   TProfile *f3pCorrelatorVsPtDiffPro = dynamic_cast<TProfile *>(listProfiles->At(5));
 233
 234   //______________________________________________________________//
 235   //Draw the histograms
 236   TCanvas *c0 = new TCanvas("c0","Analysis settings",0,0,500,500);
 237   c0->SetHighLightColor(10); c0->SetFillColor(10);
 238   fAnalysisSettings->Draw();
 239
 240   TCanvas *c1 = new TCanvas("c1","3p correlator",50,50,500,500);
 241   c1->SetHighLightColor(10); c1->SetFillColor(10);
 242   f3pCorrelatorHist->Draw("E");
 243
 244   TCanvas *c2 = new TCanvas("c2","Detector bias",100,100,500,500);
 245   c2->SetHighLightColor(10); c2->SetFillColor(10);
 246   fDetectorBiasHist->Draw("E");
 247
 248   TCanvas *c3 = new TCanvas("c3","3p correlator vs M",150,150,500,500);
 249   c3->SetHighLightColor(10); c3->SetFillColor(10);
 250   f3pCorrelatorVsMHist->Draw("E");
 251
 252   TCanvas *c4 = new TCanvas("c4","Detector bias vs M",200,200,500,500);
 253   c4->SetHighLightColor(10); c4->SetFillColor(10);
 254   fDetectorBiasVsMHist->Draw("E");
 255
 256   TCanvas *c5 = new TCanvas("c5","3p correlator (pro)",500,0,500,500);
 257   c5->SetHighLightColor(10); c5->SetFillColor(10);
 258   f3pCorrelatorPro->Draw("E");
 259
 260   TCanvas *c6 = new TCanvas("c6","Non isotropic terms",550,50,500,500);
 261   c6->SetHighLightColor(10); c6->SetFillColor(10);
 262   fNonIsotropicTermsPro->Draw("E");
 263
 264   TCanvas *c7 = new TCanvas("c7","3p correlator vs M (pro)",600,100,500,500);
 265   c7->SetHighLightColor(10); c7->SetFillColor(10);
 266   f3pCorrelatorVsMPro->Draw("E");
 267
 268   TCanvas *c8 = new TCanvas("c8","Non isotropic terms vs M (pro)",650,150,500,500);
 269   c8->SetHighLightColor(10); c8->SetFillColor(10);
 270   fNonIsotropicTermsVsMPro->Draw("colz");
 271
 272   TCanvas *c9 = new TCanvas("c9","3p correlator vs sum pt",700,200,500,500);
 273   c9->SetHighLightColor(10); c9->SetFillColor(10);
 274   f3pCorrelatorVsPtSumPro->Draw("E");
 275
 276   TCanvas *c10 = new TCanvas("c10","3p correlator vs diff pt",750,250,500,500);
 277   c10->SetHighLightColor(10); c10->SetFillColor(10);
 278   f3pCorrelatorVsPtDiffPro->Draw("E");
 279
 280   Double_t g3pCorrelatorValue = 0., g3pCorrelatorError = 0.;
 281   GetCorrelatorAndError(f3pCorrelatorVsPtSumPro,
 282                         //GetCorrelatorAndError(f3pCorrelatorVsPtDiffPro,
 283                         g3pCorrelatorValue,
 284                         g3pCorrelatorError,1,20);
 285
 286   //______________________________________________________________//
 287   //Print the results
 288   Printf("============================================================");
 289   cout<<"<cos(psi1 + psi2 - 2phi3)>: "<<
 290     g3pCorrelatorValue <<
 291     " +- " <<
 292     g3pCorrelatorError << endl;
 293   cout<<"<cos(phi1 + phi2 - 2phi3)>: "<<
 294     f3pCorrelatorHist->GetBinContent(1) <<
 295     " +- " <<
 296     f3pCorrelatorHist->GetBinError(1)<<endl;
 297   Printf("============================================================");
 298
 299   listOutput->Add(f3pCorrelatorHist);
 300   listOutput->Add(f3pCorrelatorVsMHist);
 301   listOutput->Add(f3pCorrelatorVsPtSumPro);
 302   listOutput->Add(f3pCorrelatorVsPtDiffPro);
 303
 304   return listOutput;
 305 }
 306
 307 //____________________________________________________________//
 308 void GetCorrelatorAndError(TProfile *f3pCorrelatorVsPt,
 309                            Double_t &g3pCorrelatorValue,
 310                            Double_t &g3pCorrelatorError,
 311                            Int_t iBinLow = 0,
 312                            Int_t iBinHigh = 0) {
 313   //Function to return the average value of the 3p correlator
 314   //<cos(psi1 + psi2 - 2phi3)> and its error.
 315   //The first argument is one of the 3p TProfile objects vs pt.
 316   //The second and third argument give the integral and its error.
 317   //The fourth and fifth, if specified, indicate the lowest and
 318   //highest bin the calculation should be performed for.
 319   Int_t gBinLow = 1, gBinHigh = f3pCorrelatorVsPt->GetNbinsX();
 320   if(iBinLow) gBinLow = iBinLow;
 321   if(iBinHigh) gBinHigh = iBinHigh;
 322
 323   Int_t iBinCounter = 0;
 324   Double_t gSumBinContentTimesWeight = 0., gSumWeight = 0.;
 325   Double_t gSumBinContentTimesWeightSquared = 0.;
 326   for(Int_t iBin = gBinLow; iBin <= gBinHigh; iBin++) {
 327     iBinCounter += 1;
 328
 329     gSumBinContentTimesWeight += f3pCorrelatorVsPt->GetBinContent(iBin)*f3pCorrelatorVsPt->GetBinEntries(iBin);
 330     gSumWeight += f3pCorrelatorVsPt->GetBinEntries(iBin);
 331     gSumBinContentTimesWeightSquared += TMath::Power(gSumBinContentTimesWeight,2);
 332   }
 333
 334   Printf("%lf - %d",gSumWeight,iBinCounter);
 335   //Calculate the g3pCorrelatorValue and its error
 336   g3pCorrelatorValue = -1000.;
 337   g3pCorrelatorError = 1000.;
 338   if((gSumWeight)&&(iBinCounter)) {
 339     g3pCorrelatorValue = gSumBinContentTimesWeight/(gSumWeight*iBinCounter);
 340     g3pCorrelatorError = TMath::Sqrt(gSumBinContentTimesWeightSquared/TMath::Power(gSumWeight,2))/iBinCounter;
 341   }
 342
 343   return;
 344 }
 345