PWGCF/FLOW/macros/extractJetFlow.C

   1
   2 /*
   3  * extractJetFlow.C macro
   4  *
   5  * author: Redmer Alexander Bertens, 2013
   6  * rbertens@cern.ch, r.a.bertens@uu.nl, rbertens@nikhef.nl
   7  *
   8  * this macro performs a number of calculations which are necessary to understand the output
   9  * of the PWGJE PbPb Jet Flow wagons, most notably
  10  *
  11  * [1] get delta pt widths for different background subtraction schemes
  12  * [2] extract hybrid track (jet constituent track) and jet flow
  13  *     from different flow methods
  14  */
  15
  16     // global variables
  17     TFile f("AnalysisResults.root");
  18     TFile w("extractedJetFlow.root", "RECREATE");
  19     TList* keys = f.GetListOfKeys();
  20     TList* lNoFit(0x0);         // placeholder pointer for output list no fit
  21     TList* lComb(0x0);          // placeholder pointer for output list combined fit
  22     TList* lInt(0x0);           // placeholder pointer for output list integrated v2
  23     TDirectoryFile* dirNoFit(0x0);      // placeholder dir for JetFlow
  24     TDirectoryFile* dirComb(0x0);
  25     TDirectoryFile* dirInt(0x0);
  26     // centralities
  27     Double_t _c[] = {0, 10, 20, 30, 40, 50, 60, 70, 90};
  28     TArrayD* centralities = new TArrayD((int)(sizeof(_c)/sizeof(_c[0])), _c);
  29     static const int maxCen((int)(sizeof(_c)/sizeof(_c[0])));// max number of centrality bins
  30     // pt array for jet flow analysis
  31     Double_t ptJ[] = {1, 10, 20, 30, 50, 100, 150, 200};
  32     TArrayD* _ptJ = new TArrayD(sizeof(ptJ)/sizeof(ptJ[0]), ptJ);
  33     // pt array for hybrid flow analysis
  34     Double_t ptH[] = {0., 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5};
  35     TArrayD* _ptH = new TArrayD(sizeof(ptH)/sizeof(ptH[0]), ptH);
  36     Double_t jetRadius(0.2);
  37     AliAnalysisTaskRhoVnModulation* rho = new AliAnalysisTaskRhoVnModulation();
  38     TH1F* r2V0A(0x0);      // container for second order ep resolution
  39     TH1F* r2V0C(0x0);
  40     TH1F* r3V0A(0x0);      // container for third order ep resolution
  41     TH1F* r3V0C(0x0);
  42     // placeholders for delta pt histos from LHS fit
  43     TH1F* dPtTheoryNoVn(0x0);
  44     TH1F* dPtTheoryVn(0x0);
  45     TH1F* dPtkNoFit(0x0);
  46     TH1F* dPtkComb(0x0);
  47     TH1F* dPtkInt(0x0);
  48     // placeholders for delta pt histos from RMS
  49     TH1F* RMSdPtTheoryNoVn(0x0);
  50     TH1F* RMSdPtTheoryVn(0x0);
  51     TH1F* RMSdPtkNoFit(0x0);
  52     TH1F* RMSdPtkComb(0x0);
  53     TH1F* RMSdPtkInt(0x0);
  54 //_____________________________________________________________________________
  55 void extractJetFlow() {
  56     // macro to read output of v1.0 of jet flow tasks
  57     ReadOutputFile();
  58     w.cd();
  59     // Get the delta pt info from reading RMS and Mean from the histograms
  60     if(lNoFit)  {
  61         w.mkdir(Form("DeltaPt_HISTO_%s", lNoFit->GetName()));
  62         w.cd(Form("DeltaPt_HISTO_%s", lNoFit->GetName()));
  63         RMSdPtkNoFit = GetDeltaPtRMS(lNoFit, TString("kNoFit"));
  64         RMSdPtkNoFit->Write();
  65     }
  66     if(lComb)   {
  67         w.mkdir(Form("DeltaPt_HISTO_%s", lComb->GetName()));
  68         w.cd(Form("DeltaPt_HISTO_%s", lComb->GetName()));
  69         RMSdPtkComb = GetDeltaPtRMS(lComb, TString("kComb"));
  70         RMSdPtkComb->Write();
  71     }
  72     if(lInt)    {
  73         w.mkdir(Form("DeltaPt_HISTO_%s", lInt->GetName()));
  74         w.cd(Form("DeltaPt_HISTO_%s", lInt->GetName()));
  75         RMSdPtkInt = GetDeltaPtRMS(lInt, TString("kInt"));
  76         RMSdPtkInt->Write();
  77     }
  78     // Get the delta pt info from doing a iterative LHS gaus fit
  79     if(lNoFit) {
  80         w.mkdir(Form("DeltaPt_LHSFIT_%s", lNoFit->GetName()));
  81         w.cd(Form("DeltaPt_LHSFIT_%s", lNoFit->GetName()));
  82         dPtkNoFit = GetDeltaPtSigma(lNoFit, TString("kNoFit"));
  83         dPtkNoFit->Write();
  84         GetDeltaPtMean(lNoFit, TString("kNoFit"))->Write();
  85     }
  86     if(lComb) {
  87         w.mkdir(Form("DeltaPt_LHSFIT_%s", lComb->GetName()));
  88         w.cd(Form("DeltaPt_LHSFIT_%s", lComb->GetName()));
  89         dPtkComb = GetDeltaPtSigma(lComb, TString("kComb"));
  90         dPtkComb->Write();
  91         GetDeltaPtMean(lComb, TString("kComb"))->Write();
  92     }
  93     if(lInt) {
  94         w.mkdir(Form("DeltaPt_LHSFIT_%s", lInt->GetName()));
  95         w.cd(Form("DeltaPt_LHSFIT_%s", lInt->GetName()));
  96         dPtkInt = GetDeltaPtSigma(lInt, TString("kInt"));
  97         dPtkInt->Write();
  98         GetDeltaPtMean(lInt, TString("kInt"))->Write();
  99     }
 100     // Get the delta pt predictions
 101     w.mkdir("DeltaPt_PREDICTION");
 102     GetPredictedDeltaPtSigma(lComb, "");
 103     // extract the flow
 104     TList* listNoFit = dirNoFit->GetListOfKeys();
 105     for(Int_t i(0); i < listNoFit->GetEntries(); i++) {
 106         TString string = listNoFit->At(i)->GetName();
 107         if(string.Contains("JetFlow")) {
 108             for(Int_t j(0); j < 8; j++) {
 109                 if(string.EndsWith(Form("_%i_histograms", j*10))) {
 110                     TList* op = (TList*)dirNoFit->Get(string);
 111                     GetJetTrackFlow(op, j*10-5);
 112                 }
 113             }
 114         }
 115     }
 116     // get the integrated v2 and v3 that were used for the jet
 117     // background subtraction
 118     if(lComb) GetIntegratedVn(lComb, TString("VZEROC"));
 119     // get the relative improvements
 120     GetRelativeImprovements();
 121     GetRelativeImprovementsFromRMS();
 122     // get flow of hyrbid tracks (jet constituents) and flow of jets
 123     TList* listNoFit = dirNoFit->GetListOfKeys();
 124     for(Int_t i(0); i < listNoFit->GetEntries(); i++) {
 125         TString string = listNoFit->At(i)->GetName();
 126         if(string.Contains("HybridFlow")) {
 127             for(Int_t j(0); j < 8; j++) {
 128                 if(string.EndsWith(Form("_%i_histograms", j*10))) {
 129                     TList* op = (TList*)dirNoFit->Get(string);
 130                     GetHybridTrackFlow(op, j*10-5);
 131                 }
 132             }
 133         }
 134     }
 135     TList* listComb = dirComb->GetListOfKeys();
 136     for(Int_t i(0); i < listComb->GetEntries(); i++) {
 137         TString string = listComb->At(i)->GetName();
 138         if(string.Contains("JetFlow")) {
 139             for(Int_t j(0); j < 8; j++) {
 140                 if(string.EndsWith(Form("_%i_histograms", j*10))) {
 141                     TList* op = (TList*)dirComb->Get(string);
 142                     GetJetTrackFlow(op, j*10-5);
 143                 }
 144             }
 145         }
 146     }
 147     TList* listInt = dirInt->GetListOfKeys();
 148     for(Int_t i(0); i < listInt->GetEntries(); i++) {
 149         TString string = listInt->At(i)->GetName();
 150         if(string.Contains("JetFlow")) {
 151             for(Int_t j(0); j < 8; j++) {
 152                 if(string.EndsWith(Form("_%i_histograms", j*10))) {
 153                     TList* op = (TList*)dirInt->Get(string);
 154                     GetJetTrackFlow(op, j*10-5);
 155                 }
 156             }
 157         }
 158     }
 159     // save the analysis summary histogram
 160     if(lNoFit)  GetAnalysisSummary(lNoFit);
 161     if(lComb)   GetAnalysisSummary(lComb);
 162     if(lInt)    GetAnalysisSummary(lInt);
 163     // lock and write the output file
 164     w.Close();
 165     // get started !
 166     TBrowser* browser = new TBrowser();
 167     SetStyle();
 168 }
 169 //_____________________________________________________________________________
 170 void LoadLibraries() {
 171     // Load common libraries - who knows what the future will bring ...
 172     // note that these need to be loaded prior to executing the macro (they're
 173     // necessary for the types specified as global variables)
 174     gSystem->Load("libTree");
 175     gSystem->Load("libVMC");
 176     gSystem->Load("libGeom");
 177     gSystem->Load("libGui");
 178     gSystem->Load("libXMLParser");
 179     gSystem->Load("libMinuit");
 180     gSystem->Load("libMinuit2");
 181     gSystem->Load("libProof");
 182     gSystem->Load("libPhysics");
 183     gSystem->Load("libSTEERBase");
 184     gSystem->Load("libESD");
 185     gSystem->Load("libAOD");
 186     gSystem->Load("libOADB");
 187     gSystem->Load("libANALYSIS");
 188     gSystem->Load("libCDB");
 189     gSystem->Load("libRAWDatabase");
 190   //  gSystem->Load("libSTEER");
 191     gSystem->Load("libEVGEN");
 192     gSystem->Load("libANALYSISalice");
 193     gSystem->Load("libCORRFW");
 194     gSystem->Load("libTOFbase");
 195     //gSystem->Load("libTOFrec");
 196     gSystem->Load("libRAWDatabase.so");
 197     gSystem->Load("libRAWDatarec.so");
 198     gSystem->Load("libTPCbase.so");
 199     gSystem->Load("libTPCrec.so");
 200     gSystem->Load("libITSbase.so");
 201     gSystem->Load("libITSrec.so");
 202     gSystem->Load("libTRDbase.so");
 203     gSystem->Load("libTENDER.so");
 204     gSystem->Load("libSTAT.so");
 205     gSystem->Load("libTRDrec.so");
 206     gSystem->Load("libHMPIDbase.so");
 207     gSystem->Load("libPWGPP.so");
 208     gSystem->Load("libPWGHFbase");
 209     gSystem->Load("libPWGDQdielectron");
 210     gSystem->Load("libPWGHFhfe");
 211     gSystem->Load("libEMCALUtils");
 212     gSystem->Load("libPHOSUtils");
 213     gSystem->Load("libPWGCaloTrackCorrBase");
 214     gSystem->Load("libEMCALraw");
 215     gSystem->Load("libEMCALbase");
 216     gSystem->Load("libEMCALrec");
 217     gSystem->Load("libTRDbase");
 218     gSystem->Load("libVZERObase");
 219     gSystem->Load("libVZEROrec");
 220     gSystem->Load("libTENDER");
 221     gSystem->Load("libTENDERSupplies");
 222     gSystem->Load("libPWGEMCAL");
 223     gSystem->Load("libPWGGAEMCALTasks");
 224     gSystem->Load("libPWGTools");
 225     gSystem->Load("libPWGCFCorrelationsBase");
 226     gSystem->Load("libPWGCFCorrelationsDPhi");
 227     gSystem->Load("libCGAL");
 228     gSystem->Load("libfastjet");
 229     gSystem->Load("libSISConePlugin");
 230     gSystem->Load("libJETAN");
 231     gSystem->Load("libFASTJETAN");
 232     gSystem->Load("libPWGJEEMCALJetTasks");
 233     gSystem->Load("libPWGflowBase");
 234 }
 235 //_____________________________________________________________________________
 236 void ReadOutputFile() {
 237     // loop over the keys
 238     for(Int_t i(0); i < keys->GetEntries(); i++) {
 239         if(TString(keys->At(i)->GetName()).EndsWith("PWGJE")) {
 240             if (TString(keys->At(i)->GetName()).Contains("kComb")) {
 241                 printf(" > Found PWGJE output \n \t %s \n ", keys->At(i)->GetName());
 242                 lComb = (TList*)f.Get(keys->At(i)->GetName());
 243             }
 244             if (TString(keys->At(i)->GetName()).Contains("kInt")) {
 245                 printf(" > Found PWGJE output \n \t %s \n ", keys->At(i)->GetName());
 246                 lInt = (TList*)f.Get(keys->At(i)->GetName());
 247             }
 248             if (TString(keys->At(i)->GetName()).Contains("kNoFit")) {
 249                 printf(" > Found PWGJE output \n \t %s \n ", keys->At(i)->GetName());
 250                 lNoFit = (TList*)f.Get(keys->At(i)->GetName());
 251             }
 252         }
 253         if(TString(keys->At(i)->GetName()).EndsWith("PWGCF")) {
 254             if(TString(keys->At(i)->GetName()).Contains("kComb")) {
 255                 printf(" > Found PWGCF output \n \t %s \n ", keys->At(i)->GetName());
 256                 dirComb = (TDirectoryFile*)f.Get(keys->At(i)->GetName());
 257             }
 258             if(TString(keys->At(i)->GetName()).Contains("kInt")) {
 259                 printf(" > Found PWGCF output \n \t %s \n ", keys->At(i)->GetName());
 260                 dirInt = (TDirectoryFile*)f.Get(keys->At(i)->GetName());
 261             }
 262             if(TString(keys->At(i)->GetName()).Contains("kNoFit")) {
 263                 printf(" > Found PWGCF output \n \t %s \n ", keys->At(i)->GetName());
 264                 dirNoFit = (TDirectoryFile*)f.Get(keys->At(i)->GetName());
 265             }
 266         }
 267     }   // end of loop over keys
 268 }
 269 //_____________________________________________________________________________
 270 GetRelativeImprovements() {
 271     // get the relative improvements in delta pt widths
 272     // note that the error propagation towards the relative
 273     // improvement is NOT CORRECT !
 274     if(dPtTheoryVn && dPtTheoryNoVn ) {
 275         TH1F* impTheory = new TH1F("theory, LHS fit ", "theory, LHS fit", centralities->GetSize()-1, centralities->GetArray());
 276         impTheory->GetYaxis()->SetTitle("relative improvement");
 277         impTheory->GetXaxis()->SetTitle("centrality percentile");
 278         for(Int_t i (0); i < centralities->GetSize(); i++) {
 279             Double_t a = dPtTheoryNoVn->GetBinContent(i+1);
 280             Double_t b = dPtTheoryVn->GetBinContent(i+1);
 281             if(b!=0) {
 282                 impTheory->SetBinContent(i+1, (b-a)/b);
 283                 impTheory->SetBinError(1+i, impTheory->GetBinError(i+1));
 284             }
 285         }
 286     }
 287     if(dPtkNoFit && dPtkComb) {
 288     TH1F* impComb = new TH1F("measured, LHS fit", "measured, LHS fit", centralities->GetSize()-1, centralities->GetArray());
 289         impComb->GetYaxis()->SetTitle("relative improvement");
 290         impComb->GetXaxis()->SetTitle("centrality percentile");
 291         for(Int_t i (0); i < centralities->GetSize(); i++) {
 292             Double_t a = dPtkComb->GetBinContent(i+1);
 293             Double_t b = dPtkNoFit->GetBinContent(i+1);
 294             if(b!=0) {
 295                 impComb->SetBinContent(i+1, (b-a)/b);
 296                 impComb->SetBinError(i+1, impComb->GetBinError(i+1));
 297             }
 298         }
 299     }
 300     if(dPtkNoFit && dPtkInt) {
 301     TH1F* impInt = new TH1F("relative improvement #delta p_{T} #sigma, kInt ", "relative improvement #delta p_{T} #sigma, kInt", centralities->GetSize()-1, centralities->GetArray());
 302         impInt->GetYaxis()->SetTitle("relative improvement");
 303         impInt->GetXaxis()->SetTitle("centrality percentile");
 304         for(Int_t i (0); i < centralities->GetSize(); i++) {
 305             Double_t a = dPtkInt->GetBinContent(i+1);
 306             Double_t b = dPtkNoFit->GetBinContent(i+1);
 307             if(b!=0) {
 308                 impInt->SetBinContent(i+1, (b-a)/b);
 309                 impInt->SetBinError(1+i, impInt->GetBinError(i+1));
 310             }
 311         }
 312     }
 313     // write the output to file
 314     w.cd();
 315     w.mkdir("Relative improvement delta pt distributions");
 316     w.cd("Relative improvement delta pt distributions");
 317     FormatMe(impTheory);
 318     impTheory->Write();
 319     FormatMe(impComb);
 320     impComb->Write();
 321     FormatMe(impInt);
 322     impInt->Write();
 323
 324 }
 325 //_____________________________________________________________________________
 326 GetRelativeImprovementsFromRMS() {
 327     // get the relative improvements in delta pt widths
 328     // note that the error propagation towards the relative
 329     // improvement is NOT CORRECT !
 330     if(dPtTheoryVn && dPtTheoryNoVn ) {
 331         TH1F* impTheory = new TH1F("theory", "theory", centralities->GetSize()-1, centralities->GetArray());
 332         impTheory->GetYaxis()->SetTitle("relative improvement");
 333         impTheory->GetXaxis()->SetTitle("centrality percentile");
 334         for(Int_t i (0); i < centralities->GetSize(); i++) {
 335             Double_t a = RMSdPtTheoryNoVn->GetBinContent(i+1);
 336             Double_t b = RMSdPtTheoryVn->GetBinContent(i+1);
 337             if(b!=0) {
 338                 impTheory->SetBinContent(i+1, (b-a)/b);
 339                 impTheory->SetBinError(1+i, impTheory->GetBinError(i+1));
 340             }
 341         }
 342     }
 343     if(RMSdPtkNoFit && RMSdPtkComb) {
 344     TH1F* impComb = new TH1F("measured", "measured", centralities->GetSize()-1, centralities->GetArray());
 345         impComb->GetYaxis()->SetTitle("relative improvement");
 346         impComb->GetXaxis()->SetTitle("centrality percentile");
 347         for(Int_t i (0); i < centralities->GetSize(); i++) {
 348             Double_t a = RMSdPtkComb->GetBinContent(i+1);
 349             Double_t b = RMSdPtkNoFit->GetBinContent(i+1);
 350             if(b!=0) {
 351                 impComb->SetBinContent(i+1, (b-a)/b);
 352                 impComb->SetBinError(i+1, impComb->GetBinError(i+1));
 353             }
 354         }
 355     }
 356     if(RMSdPtkNoFit && RMSdPtkInt) {
 357     TH1F* impInt = new TH1F("measured ", "measured ", centralities->GetSize()-1, centralities->GetArray());
 358         impInt->GetYaxis()->SetTitle("relative improvement");
 359         impInt->GetXaxis()->SetTitle("centrality percentile");
 360         for(Int_t i (0); i < centralities->GetSize(); i++) {
 361             Double_t a = RMSdPtkInt->GetBinContent(i+1);
 362             Double_t b = RMSdPtkNoFit->GetBinContent(i+1);
 363             if(b!=0) {
 364                 impInt->SetBinContent(i+1, (b-a)/b);
 365                 impInt->SetBinError(1+i, impInt->GetBinError(i+1));
 366             }
 367         }
 368     }
 369     // write the output to file
 370     w.cd();
 371     w.mkdir("Relative improvement delta pt distributions from RMS");
 372     w.cd("Relative improvement delta pt distributions from RMS");
 373     FormatMe(impTheory);
 374     impTheory->Write();
 375     FormatMe(impComb);
 376     impComb->Write();
 377     FormatMe(impInt);
 378     impInt->Write();
 379
 380 }
 381 //_____________________________________________________________________________
 382 void GetHybridTrackFlow(TList* jf, Int_t c) {
 383     // get hybrid track flow from the vzero ep and qc anlaysis
 384     if(!jf) {
 385         printf(" > couldn't find output list with name %s \n", name.Data());
 386         return;
 387     }
 388     w.mkdir(Form("GetHybridTrackFlow_%s", jf->GetName()));
 389     w.cd(Form("GetHybridTrackFlow_%s", jf->GetName()));
 390     TProfile* v0a = (TProfile*)jf->FindObject("Differential v_{2}^{obs} VZEROA");
 391     TProfile* v0c = (TProfile*)jf->FindObject("Differential v_{2}^{obs} VZEROC");
 392     TProfile* qc2 = (TProfile*)jf->FindObject("Differential cumulants v2");
 393     TProfile* rc  = (TProfile*)jf->FindObject("Reference cumulants");
 394     if(qc2 && rc && v0a ) {
 395         TH1F* result = rho->GetDifferentialQC(rc, qc2, _ptH, 2);
 396         FormatMe(result);
 397         TString t = "qc2_";
 398         t+=jf->GetName();
 399         result->SetNameTitle(t.Data(), t.Data());
 400         for(Int_t i(0); i < result->GetXaxis()->GetNbins(); i++) {
 401 //            Double_t M(rc->GetBinEntries(1)/2./700.);
 402 //            Double_t Mp(qc2->GetBinEntries(1+i)/qc2->GetEntries());
 403 //            Double_t errinv(rc->GetBinContent(1)*TMath::Sqrt(M*Mp));
 404 //            cout << " errinv " << errinv << endl;
 405 //            if(errinv > 0) errinv = TMath::Sqrt(errinv);
 406 //            cout << " err " << errinv << endl;
 407             result->SetBinError(1+i, v0a->GetBinError(1+i));
 408         }
 409         result->Write();
 410     }
 411     if(v0a) {
 412         TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0a, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, c, 2);
 413         FormatMe(result);
 414         result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 415         result->GetYaxis()->SetTitle("v_{2}");
 416         result->Write();
 417     }
 418     if(v0c) {
 419         TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0c, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, c, 2);
 420         FormatMe(result);
 421         result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 422         result->GetYaxis()->SetTitle("v_{2}");
 423        result->Write();
 424     }
 425     // attempt to get the flow from the qc analysis
 426     TDirectoryFile* qc = (TDirectoryFile*)f.Get("QC");
 427     if(qc) {
 428         TList* qcl = (TList*)qc->Get(Form("QC_hybrid_flow_%s", app.Data()));
 429         if(qcl) {
 430             AliFlowCommonHistResults* flow = (AliFlowCommonHistResults*)qcl->FindObject("AliFlowCommonHistResults2ndOrderQC");
 431             if(flow) flow->GetHistDiffFlowPtPOI()->Write();
 432         }
 433     }
 434 }
 435 //_____________________________________________________________________________
 436 void GetJetTrackFlow(TList* jf, Int_t c) {
 437     // get hybrid track flow from the vzero ep and qc anlaysis
 438     if(!jf) {
 439         printf(" > couldn't find output list with name %s \n", name.Data());
 440         return;
 441     }
 442     w.mkdir(Form("GetJetTrackFlow_%s", jf->GetName()));
 443     w.cd(Form("GetJetTrackFlow_%s", jf->GetName()));
 444     TProfile* v0a = (TProfile*)jf->FindObject("Differential v_{2}^{obs} VZEROA");
 445     TProfile* v0c = (TProfile*)jf->FindObject("Differential v_{2}^{obs} VZEROC");
 446     TProfile* qc2 = (TProfile*)jf->FindObject("Differential cumulants v2");
 447     TProfile* rc  = (TProfile*)jf->FindObject("Reference cumulants");
 448         if(qc2 && rc ) {
 449         TH1F* result = rho->GetDifferentialQC(rc, qc2, _ptJ, 2);
 450         TString t = "qc2_";
 451         t+=jf->GetName();
 452          result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 453         result->GetYaxis()->SetTitle("v_{2}");
 454        result->SetNameTitle(t.Data(), t.Data());
 455         FormatMe(result);
 456         result->Write();
 457     }
 458     if(v0a) {
 459         TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0a, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, c, 2);
 460          result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 461         result->GetYaxis()->SetTitle("v_{2}");
 462        FormatMe(result);
 463           result->Write();
 464     }
 465     if(v0c) {
 466         TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0c, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, c, 2);
 467          result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 468         result->GetYaxis()->SetTitle("v_{2}");
 469        FormatMe(result);
 470         result->Write();
 471     }
 472 }
 473 //_____________________________________________________________________________
 474 TH1F* GetDeltaPtRMS(TList* l, TString suffix) {
 475     // get the RMS value of delta pt
 476     TH1F* deltaPtRMS = new TH1F(Form("#delta p_{T} RMS, %s", suffix.Data()), Form("#delta p_{T} RMS, %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
 477     deltaPtRMS->GetXaxis()->SetTitle("centrality percentile");
 478     deltaPtRMS->GetYaxis()->SetTitle("RMS [GeV/c]");
 479     for(Int_t i(0); i < maxCen; i++) {
 480         TString name = Form("fHistDeltaPtDeltaPhi2_%i", i);
 481         printf(" > searching for %s in list \n %s < \n ", name.Data(), l->GetName());
 482         TH2D* dpt((TH2D*)l->FindObject(name.Data()));
 483         if(!dpt) continue;
 484         deltaPtRMS->SetBinContent(i+1, dpt->GetRMS(2));
 485         deltaPtRMS->SetBinError(i+1, dpt->GetRMSError(2));
 486     }
 487     FormatMe(deltaPtRMS);
 488     return deltaPtRMS;
 489 }
 490 //_____________________________________________________________________________
 491 TH1F* GetDeltaPtSigma(TList* l, TString suffix) {
 492     // get the sigma of the delta pt distribution from a recursive LHS gauss fit
 493     TH1F* deltaPtMean = new TH1F(Form("#delta p_{T} #sigma, %s", suffix.Data()), Form("#delta p_{T} #sigma, %s",suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
 494     deltaPtMean->GetYaxis()->SetTitle("#sigma [GeV/c]");
 495     deltaPtMean->GetXaxis()->SetTitle("centrality percentile");
 496     for(Int_t i(0); i < maxCen; i++) {
 497         TString name = Form("fHistDeltaPtDeltaPhi2_%i", i);
 498         printf(" > searching for %s in list \n \t  %s < \n ", name.Data(), l->GetName());
 499         TH2D* dpt((TH2D*)l->FindObject(name.Data()));
 500         if(!dpt) continue;
 501         TH1D* temp = dpt->ProjectionY(/*"_py", 0, -1, "e"*/);    // do error propagation for projection
 502         Double_t s = temp->GetRMS();
 503         Double_t m = temp->GetMean();
 504         TF1* fit = new TF1(Form("sigma_%s", temp->GetName()), "gaus", m-3*s, m+0.5*s);
 505         for(Int_t j(0); j < 10; j++) {
 506             Double_t _m(m), _s(s);
 507             temp->Fit(fit, "QILR");
 508             fit->SetRange(m-3*s, m+0.5*s);
 509             m = fit->GetParameter(1);
 510             s = fit->GetParameter(2);
 511
 512         }
 513         deltaPtMean->SetBinContent(1+i, fit->GetParameter(2));
 514         deltaPtMean->SetBinError(1+i, fit->GetParError(2));
 515     }
 516     FormatMe(deltaPtMean);
 517     return deltaPtMean;
 518 }
 519 //_____________________________________________________________________________
 520 TH1F* GetDeltaPtMean(TList* l, TString suffix) {
 521     // get the mean of the delta pt distribution from a recursive LHS gauss fit
 522     TH1F* deltaPtMean = new TH1F(Form("#delta p_{T} mean %s", suffix.Data()), Form("#delta p_{T} mean %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
 523     deltaPtMean->GetYaxis()->SetTitle("mean [GeV/c]");
 524     deltaPtMean->GetXaxis()->SetTitle("centrality percentile");
 525     for(Int_t i(0); i < maxCen; i++) {
 526         TString name = Form("fHistDeltaPtDeltaPhi2_%i", i);
 527         printf(" > searching for %s in list \n \t %s < \n ", name.Data(), l->GetName());
 528         TH2D* dpt((TH2D*)l->FindObject(name.Data()));
 529         if(!dpt) continue;
 530         TH1D* temp = dpt->ProjectionY(/*"_py", 0, -1, "e"*/);    // do error propagation for projection
 531         Double_t s = temp->GetRMS();
 532         Double_t m = temp->GetMean();
 533         TF1* fit = new TF1(Form("mean_%s", temp->GetName()), "gaus", m-3*s, m+0.5*s);
 534         TH1F* qam = new TH1F(Form("mu_%s", temp->GetName()), "#mu_{i} / #mu_{i-1}", 10, 0, 10);
 535         TH1F* qas = new TH1F(Form("sigma_%s", temp->GetName()), "#sigma_{i} / #sigma_{i-1}", 10, 0, 10);
 536         fit->SetParLimits(2, s/2., s*2.);
 537         for(Int_t j(0); j < 10; j++) {
 538             Double_t _m(m), _s(s);
 539             fit->SetRange(m-3*s, m+0.5*s);
 540             temp->Fit(fit, "QILR");
 541             m = fit->GetParameter(1);
 542             s = fit->GetParameter(2);
 543             if(!m == 0) qam->SetBinContent(j+1, _m/m);
 544             if(!s == 0) qas->SetBinContent(j+1, _s/s);
 545         }
 546         deltaPtMean->SetBinContent(1+i, fit->GetParameter(1));
 547         deltaPtMean->SetBinError(1+i, fit->GetParError(1));
 548         temp->Write();
 549         qas->Write();
 550         qam->Write();
 551     }
 552     FormatMe(deltaPtMean);
 553     return deltaPtMean;
 554 }
 555 //_____________________________________________________________________________
 556 void GetPredictedDeltaPtSigma(TList* l, TString suffix) {
 557     // get predicted delta pt sigma
 558     TH1F* deltaPtSigma = new TH1F(Form("predicted #delta p_{T} #sigma %s", suffix.Data()), Form("predicted #delta p_{T} #sigma %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
 559     deltaPtSigma->GetYaxis()->SetTitle("predicted #sigma [GeV/c]");
 560     deltaPtSigma->GetXaxis()->SetTitle("centrality percentile");
 561     TH1F* deltaPtSigmaNoV = new TH1F("predicted #delta p_{T} #sigma no vn", "predicted #delta p_{T} #sigma no vn", centralities->GetSize()-1, centralities->GetArray());
 562     deltaPtSigmaNoV->GetYaxis()->SetTitle("predicted #sigma [GeV/c]");
 563     deltaPtSigmaNoV->GetXaxis()->SetTitle("centrality percentile");
 564     // get the info from the methods of the class
 565     rho->SetOutputList((TList*)l->Clone());
 566     // get the resolution for the desired detector
 567     r2V0A = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, 2, centralities);
 568     r2V0A->SetNameTitle("VZEROA resolution for #Psi_{2}", "VZEROA resolution for #Psi_{2}");
 569     r3V0A = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, 3, centralities);
 570     r3V0A->SetNameTitle("VZEROA resoltuion for #Psi_{3}", "VZEROA resolution for #Psi_{3}");
 571     r2V0C = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, 2, centralities);
 572     r2V0C->SetNameTitle("VZEROC resolution for #Psi_{2}", "VZEROC resolution for #Psi_{2}");
 573     r3V0C = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, 3, centralities);
 574     r3V0C->SetNameTitle("VZEROC resolution for #Psi_{3}", "VZEROC resolution for #Psi_{3}");
 575     // grab the v2 and v3 values and do a resolution correction
 576     TH1F* v2 = new TH1F("v2", "v2", centralities->GetSize()-1, centralities->GetArray());
 577     TH1F* v3 = new TH1F("v3", "v3", centralities->GetSize()-1, centralities->GetArray());
 578     TProfile* pv2 = (TProfile*)l->FindObject("fProfV2");
 579     TProfile* pv3 = (TProfile*)l->FindObject("fProfV3");
 580     for(Int_t i(0); i < maxCen-1; i++) {
 581         v2->SetBinContent(1+i, pv2->GetBinContent(1+i));
 582         v2->SetBinError(1+i, pv2->GetBinError(1+i));
 583         v3->SetBinContent(1+i, pv3->GetBinContent(1+i));
 584         v3->SetBinError(1+i, pv3->GetBinError(1+i));
 585     }
 586     // from
 587     TH1F* cv2 = new TH1F("v2 in from literaturet","v2 int from literature",10,0,100);
 588     Double_t c_v2[] = {0, 0.03565825, 0.06394614, 0.08306863, 0.09470311, 0.09927855, 0.09630484, 0.08708335,   0.07051519, 0, 0};
 589     cv2->SetContent(c_v2);
 590     TH1F* cv3 = new TH1F("v3 int from literature","v3 int from literature", 10, 0, 100);
 591     Double_t c_v3[] = {0, 0.02159083, 0.02642751, 0.02928424, 0.03052121, 0.03004316, 0, 0, 0, 0, 0};
 592     cv3->SetContent(c_v3);
 593     for(Int_t i(0); i < centralities->GetSize()-1; i++) {
 594         TH1F* h = (TH1F*)l->FindObject(Form("fHistPicoTrackMult_%i", i));
 595         TH1F* m = (TH1F*)l->FindObject(Form("fHistPicoTrackPt_%i", i));
 596         if(!h||!m) {
 597             printf(" > Panic, one of the pt histos not found ! < \n");
 598             continue;
 599         }
 600         Double_t na(h->GetMean()*jetRadius*jetRadius*TMath::Pi()/(2.*.9*TMath::TwoPi()));
 601         Double_t naErr(h->GetMeanError()*jetRadius*jetRadius*TMath::Pi()/(2.*.9*TMath::TwoPi()));
 602         Double_t totErr(TMath::Sqrt((m->GetRMS()*m->GetRMS()+m->GetMean()*m->GetMean())*(m->GetRMS()*m->GetRMS()+m->GetMean()*m->GetMean())*naErr*naErr+4.*na*na*m->GetRMS()*m->GetRMS()*m->GetRMSError()*m->GetRMSError()+4.*na*na*m->GetMean()*m->GetMean()*m->GetMeanError()*m->GetMeanError()));
 603         Double_t err(m->GetRMS());
 604         Double_t mean(m->GetMean());
 605         Double_t _v2(cv2->GetBinContent(1+i));
 606         Double_t _v3(cv3->GetBinContent(1+i));
 607         // no v2
 608         Double_t dptnovn = TMath::Sqrt(na*(err*err+mean*mean));
 609         Double_t dpt = TMath::Sqrt(na*err*err+(na+2.*na*na*(_v2*_v2+_v3*_v3))*mean*mean);
 610         deltaPtSigma->SetBinContent(1+i, dpt);
 611         deltaPtSigma->SetBinError(1+i, totErr);
 612         deltaPtSigmaNoV->SetBinContent(1+i, dptnovn);
 613         deltaPtSigmaNoV->SetBinError(1+i, totErr);       // estimate without vn
 614     }
 615     w.mkdir("RhoTaskVnEstimates");
 616     w.cd("RhoTaskVnEstimates");
 617     FormatMe(r2V0A);
 618     r2V0A->Write();
 619     FormatMe(r3V0A);
 620     r3V0A->Write();
 621     FormatMe(r2V0C);
 622     r2V0C->Write();
 623     FormatMe(r3V0C);
 624     r3V0C->Write();
 625     FormatMe(cv2);
 626     cv2->Write();
 627     FormatMe(cv3);
 628     cv3->Write();
 629     w.cd("DeltaPt_PREDICTION");
 630     dPtTheoryVn = deltaPtSigma;
 631     RMSdPtTheoryVn = deltaPtSigma;
 632     FormatMe(dPtTheoryVn);
 633     dPtTheoryVn->Write();
 634     dPtTheoryNoVn = deltaPtSigmaNoV;
 635     RMSdPtTheoryNoVn = deltaPtSigmaNoV;
 636     FormatMe(dPtTheoryNoVn);
 637     dPtTheoryNoVn->Write();
 638 }
 639 //_____________________________________________________________________________
 640 void GetIntegratedVn(TList* l, TString det = "VZEROC") {
 641     // get the v2 and v3 values that were used to estimate local energy denstity
 642     w.cd("RhoTaskVnEstimates");
 643     TH1F* v2 = new TH1F("v2obs", "v2obs", centralities->GetSize()-1, centralities->GetArray());
 644     TH1F* v3 = new TH1F("v3obs", "v3obs", centralities->GetSize()-1, centralities->GetArray());
 645     TH1F* cv2(0x0);
 646     TH1F* cv3(0x0);
 647     TProfile* pv2 = (TProfile*)l->FindObject("fProfV2");
 648     TProfile* pv3 = (TProfile*)l->FindObject("fProfV3");
 649     for(Int_t i(0); i < maxCen; i++) {
 650         v2->SetBinContent(1+i, pv2->GetBinContent(1+i));
 651         v2->SetBinError(1+i, pv2->GetBinError(1+i));
 652         v3->SetBinContent(1+i, pv3->GetBinContent(1+i));
 653         v3->SetBinError(1+i, pv3->GetBinError(1+i));
 654     }
 655     if(det.EqualTo("VZEROA")) {
 656         cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, 2);
 657         cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, 3);
 658     } else if (det.EqualTo("VZEROC")) {
 659         cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, 2);
 660         cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, 3);
 661     } else if (det.EqualTo("VZEROComb")) {
 662         cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROComb, centralities, 2);
 663         cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROComb, centralities, 3);
 664     }
 665     TString nt2 = Form("v2, %s", det.Data());
 666     TString nt3 = Form("v3, %s", det.Data());
 667     cv2->SetNameTitle(nt2.Data(), nt2.Data());
 668     cv3->SetNameTitle(nt3.Data(), nt3.Data());
 669     FormatMe(cv2);
 670     cv2->Write();
 671     FormatMe(cv3);
 672     cv3->Write();
 673 }
 674 //_____________________________________________________________________________
 675 void GetAnalysisSummary(TList* l) {
 676     // get and format the analyis summary histogram
 677     TH1F* h = (TH1F*)l->FindObject("fHistAnalysisSummary");
 678     if(h) {
 679         Double_t iter = h->GetBinContent(37);
 680         if(iter <= 0) return;   // zero events in sample ...
 681         Int_t type = TMath::Nint(h->GetBinContent(34)/iter);
 682         TString  name = "";
 683         if(type==0) name+="kNoFit";
 684         if(type==1) name+="kV2";
 685         if(type==2) name+="kV3";
 686         if(type==3) name+="kCombined";
 687         if(type==4) name+="kFourierSeries";
 688         if(type==5) name+="kIntegratedFlow";
 689         if(type==6) name+="kQC2";
 690         if(type==7) name+="kQC4";
 691         for(Int_t i(0); i < h->GetXaxis()->GetNbins(); i++) h->SetBinContent(i+1, h->GetBinContent(i+1)/iter);
 692         w.mkdir(Form("Summary_%s", name.Data()));
 693         w.cd(Form("Summary_%s", name.Data()));
 694         h->Write();
 695     }
 696 }
 697 //_____________________________________________________________________________
 698 void SetStyle() {
 699     // set global style
 700     gStyle->SetOptStat(0);
 701     gPad->SetGrid(1,1);
 702     gPad->SetTicks(1,1);
 703     gStyle->ToggleEditor();
 704 }
 705 //_____________________________________________________________________________
 706 TH1* FormatMe(TObject* object, Int_t color = -1) {
 707     if(color=-1) color = TMath::Nint(gRandom->Uniform(1, 9));
 708     TH1* dud = (dynamic_cast<TH1*>object);
 709     if (dud) return FormatHistogram(dud, color);
 710 }
 711 //_____________________________________________________________________________
 712 TH1F* FormatHistogram(TH1* hist, Int_t color) {
 713     // return a more readable TH1F
 714     hist->SetLineWidth(3);
 715     hist->SetLineColor(color);
 716     hist->SetMarkerStyle(20);
 717     hist->SetMarkerColor(color);
 718     hist->SetMarkerColor(color);
 719     TString name = Form("%s R = %.2f", hist->GetTitle(), jetRadius);
 720     hist->SetNameTitle(name.Data(), name.Data());
 721
 722 }
 723 //_____________________________________________________________________________