PWG/muon/AliAnalysisTaskFlowSingleMu.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliAnalysisTaskFlowSingleMu.cxx 55545 2012-04-04 07:16:39Z pcrochet $ */
  17
  18 //-----------------------------------------------------------------------------
  19 /// \class AliAnalysisTaskFlowSingleMu
  20 /// Analysis task for flow of single muons in the spectrometer.
  21 /// The output is a list of histograms and CF containers.
  22 /// The macro class can run on AODs or ESDs.
  23 /// If Monte Carlo information is present, some basics checks are performed.
  24 ///
  25 /// \author Diego Stocco
  26 //-----------------------------------------------------------------------------
  27
  28 #define AliAnalysisTaskFlowSingleMu_cxx
  29
  30 #include "AliAnalysisTaskFlowSingleMu.h"
  31
  32 // ROOT includes
  33 #include "TROOT.h"
  34 #include "TH1.h"
  35 #include "TH2.h"
  36 #include "TGraphErrors.h"
  37 #include "TAxis.h"
  38 #include "TCanvas.h"
  39 #include "TLegend.h"
  40 #include "TMath.h"
  41 #include "TObjString.h"
  42 #include "TObjArray.h"
  43 #include "TF1.h"
  44 #include "TStyle.h"
  45 //#include "TMCProcess.h"
  46 #include "TFitResultPtr.h"
  47 #include "TFile.h"
  48 #include "TArrayI.h"
  49 #include "TArrayD.h"
  50 #include "TSystem.h"
  51 #include "TGraphErrors.h"
  52 #include "TLatex.h"
  53
  54 // STEER includes
  55 #include "AliAODEvent.h"
  56 #include "AliAODTrack.h"
  57 #include "AliAODMCParticle.h"
  58 #include "AliMCEvent.h"
  59 #include "AliMCParticle.h"
  60 #include "AliESDEvent.h"
  61 #include "AliESDMuonTrack.h"
  62 #include "AliVHeader.h"
  63 #include "AliAODMCHeader.h"
  64 #include "AliStack.h"
  65 #include "AliEventplane.h"
  66
  67 // ANALYSIS includes
  68 #include "AliAnalysisManager.h"
  69
  70 // CORRFW includes
  71 #include "AliCFContainer.h"
  72 #include "AliCFGridSparse.h"
  73 #include "AliCFEffGrid.h"
  74
  75 // PWG includes
  76 #include "AliVAnalysisMuon.h"
  77 #include "AliMergeableCollection.h"
  78 #include "AliCounterCollection.h"
  79 #include "AliMuonEventCuts.h"
  80 #include "AliMuonTrackCuts.h"
  81 #include "AliAnalysisMuonUtility.h"
  82
  83
  84 /// \cond CLASSIMP
  85 ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
  86 /// \endcond
  87
  88 using std::ifstream;
  89
  90 //________________________________________________________________________
  91 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
  92   AliVAnalysisMuon(),
  93   fEPKeys(0x0),
  94   fRandom(0x0)
  95 {
  96   /// Default ctor.
  97 }
  98
  99 //________________________________________________________________________
 100 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
 101   AliVAnalysisMuon(name, cuts),
 102   fEPKeys(0x0),
 103   fRandom(0x0)
 104 {
 105   //
 106   /// Constructor.
 107   //
 108   TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
 109   fEPKeys = epMethods.Tokenize(" ");
 110 }
 111
 112
 113 //________________________________________________________________________
 114 AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
 115 {
 116   //
 117   /// Destructor
 118   //
 119   delete fEPKeys;
 120   delete fRandom;
 121 }
 122
 123 //___________________________________________________________________________
 124 void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
 125 {
 126   //
 127   /// Create outputs
 128   //
 129
 130   if ( ! fRandom ) fRandom = new TRandom3();
 131   // Set seed here, to avoid having the same seed for all jobs in grid
 132   fRandom->SetSeed(0);
 133
 134   Int_t nPtBins = 160;
 135   Double_t ptMin = 0., ptMax = 80.;
 136   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
 137
 138   Int_t nEtaBins = 25;
 139   Double_t etaMin = -4.5, etaMax = -2.;
 140   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
 141
 142   Int_t nPhiBins = 36;
 143   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
 144   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
 145
 146   Int_t nDeltaPhiBins = 18;
 147   Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
 148   TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
 149
 150   Int_t nChargeBins = 2;
 151   Double_t chargeMin = -2, chargeMax = 2.;
 152   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
 153
 154   Int_t nMotherTypeBins = kNtrackSources;
 155   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
 156   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
 157
 158   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
 159   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
 160   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
 161   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
 162   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
 163
 164   Int_t iobj=0;
 165   TH1* histo = 0x0;
 166   TString currTitle = "";
 167   for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
 168     TString currEP = fEPKeys->At(iep)->GetName();
 169     AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
 170
 171     TString currName = "";
 172     for ( Int_t idim = 0; idim<kNvars; idim++){
 173       currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
 174       currName.ReplaceAll("()","");
 175
 176       cfContainer->SetVarTitle(idim, currName.Data());
 177       cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
 178     }
 179
 180     TString stepTitle[kNsteps] = {"reconstructed", "generated"};
 181
 182     TAxis* currAxis = 0x0;
 183     for (Int_t istep=0; istep<kNsteps; istep++){
 184       cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
 185       AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
 186
 187       currAxis = gridSparse->GetAxis(kHvarMotherType);
 188       for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
 189         currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
 190       }
 191     }
 192
 193     AddObjectToCollection(cfContainer, iobj++);
 194
 195     Bool_t isTPC = currEP.Contains("TPC");
 196     Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
 197     Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
 198
 199     histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
 200     histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
 201     AddObjectToCollection(histo, iobj++);
 202
 203     for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
 204       TString auxEP = fEPKeys->At(jep)->GetName();
 205       currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
 206       histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
 207       histo->GetXaxis()->SetTitle(currTitle.Data());
 208       AddObjectToCollection(histo, iobj++);
 209     }
 210
 211     currTitle = Form("cos(#psi_{1} - #psi_{2})");
 212     histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
 213     histo->GetXaxis()->SetTitle(currTitle.Data());
 214
 215     AddObjectToCollection(histo, iobj++);
 216   } // loop on EPs
 217
 218   fMuonTrackCuts->Print("mask");
 219 }
 220
 221 //________________________________________________________________________
 222 void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
 223 {
 224   //
 225   /// Fill output objects
 226   //
 227
 228   //
 229   // Base event cuts
 230   //
 231   //Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes
 232   //if ( isPileupFromSPD ) return;
 233
 234   TArrayD psiPlane(fEPKeys->GetEntriesFast());
 235   psiPlane.Reset(-999.);
 236   const Double_t kLowestVal = -10.;
 237   Int_t harmonic = 2;
 238
 239   //
 240   // Fill EP info
 241   //
 242   AliEventplane* eventPlane = InputEvent()->GetEventplane();
 243   psiPlane[0] = eventPlane->GetEventplane("V0A",InputEvent(),harmonic);
 244   psiPlane[1] = eventPlane->GetEventplane("V0C",InputEvent(),harmonic);
 245   psiPlane[2] = eventPlane->GetEventplane("Q");
 246   if ( psiPlane[2] < 0.) psiPlane[2] = -999.;
 247   Double_t qx = 0., qy = 0.;
 248   psiPlane[3] = eventPlane->CalculateVZEROEventPlane(InputEvent(),0,harmonic,qx,qy); // V0C ring 1
 249   psiPlane[4] = eventPlane->CalculateVZEROEventPlane(InputEvent(),3,harmonic,qx,qy); // V0C ring 3
 250   //psiPlane[3] = fRandom->Uniform(-TMath::Pi()/2., TMath::Pi()/2.);
 251
 252   Bool_t noValidEP = kTRUE;
 253   TArrayD psiPlaneCalc(fEPKeys->GetEntriesFast());
 254   for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
 255     psiPlaneCalc[iep] = psiPlane[iep];
 256     if ( psiPlane[iep] < kLowestVal ) continue;
 257     if ( psiPlane[iep] < 0. ) psiPlaneCalc[iep] += TMath::Pi();
 258     noValidEP = kFALSE;
 259   }
 260   if ( noValidEP ) return;
 261
 262   for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
 263     if ( psiPlane[iep] < kLowestVal ) continue;
 264     for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 265       TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 266       TString currEP = fEPKeys->At(iep)->GetName();
 267       ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hEventPlane%s",currEP.Data())))->Fill(psiPlane[iep]);
 268       for ( Int_t jep=iep+1; jep<fEPKeys->GetEntriesFast(); jep++ ) {
 269         if ( psiPlane[jep] < kLowestVal ) continue;
 270         ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub3%s%s",currEP.Data(), fEPKeys->At(jep)->GetName())))->Fill(TMath::Cos(2.*(psiPlaneCalc[iep]-psiPlaneCalc[jep])));
 271       } // loop on auxialiary EP
 272       if ( iep == 2 ) {
 273         TVector2* qsub1 = eventPlane->GetQsub1();
 274         TVector2* qsub2 = eventPlane->GetQsub2();
 275         if ( qsub1 && qsub2 )
 276           ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.));
 277       }
 278       //else ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(2.*psiPlane[iep])); // REMEMBER TO CUT
 279     } // loop on trigger
 280   } // loop on EP
 281
 282   //
 283   // Fill track info
 284   //
 285   Double_t containerInput[kNvars];
 286   AliVParticle* track = 0x0;
 287   Int_t nSteps = MCEvent() ? 2 : 1;
 288   for ( Int_t istep = 0; istep<nSteps; ++istep ) {
 289     Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks();
 290     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
 291       track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
 292
 293       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
 294       if ( ! isSelected ) continue;
 295
 296       Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
 297
 298
 299       containerInput[kHvarPt]         = track->Pt();
 300       containerInput[kHvarEta]        = track->Eta();
 301       containerInput[kHvarPhi]        = track->Phi();
 302       containerInput[kHvarCharge]     = track->Charge()/3.;
 303       containerInput[kHvarMotherType] = (Double_t)trackSrc;
 304
 305       for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
 306         if ( psiPlane[iep] < kLowestVal ) continue;
 307         TString currEP = fEPKeys->At(iep)->GetName();
 308         Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
 309
 310         // The difference between phi - psiPlane must be between (0, pi)
 311         // Since the reaction plain is symmetric in pi,
 312         // we can do in such a way that:
 313         // -pi < delta_phi + N*pi < pi
 314         // We choose N accrodingly
 315         Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
 316         deltaPhi += nPi * TMath::Pi();
 317         containerInput[kHvarDeltaPhi] = deltaPhi;
 318
 319         for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 320           TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 321           if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
 322           ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
 323         } // loop on selected trigger classes
 324       } // loop on event plane
 325     } // loop on tracks
 326   } // loop on container steps
 327 }
 328
 329
 330 //________________________________________________________________________
 331 TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
 332 {
 333   //
 334   /// Get the range from string
 335   //
 336   TArrayD centralityRange(2);
 337   centralityRange.Reset(-999.);
 338   TObjArray* array = sRange.Tokenize("_");
 339   if ( array->GetEntries() == 2 ) {
 340     for ( Int_t ientry=0; ientry<2; ientry++ ) {
 341       centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
 342     }
 343   }
 344   delete array;
 345   return centralityRange;
 346 }
 347
 348 //________________________________________________________________________
 349 void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
 350   //
 351   /// Draw some histograms at the end.
 352   //
 353
 354   AliVAnalysisMuon::Terminate("");
 355
 356   if ( ! fMergeableCollection ) return;
 357
 358   TString physSel = fTerminateOptions->At(0)->GetName();
 359   TString trigClassName = fTerminateOptions->At(1)->GetName();
 360   TString centralityRange = fTerminateOptions->At(2)->GetName();
 361   TString furtherOpt = fTerminateOptions->At(3)->GetName();
 362
 363   TObjArray* centralityClasses = centralityRange.Tokenize(",");
 364   TArrayD centralityArray(centralityClasses->GetEntries()+1);
 365   Int_t ib = 0;
 366   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 367     TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
 368     if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
 369     centralityArray.AddAt(range[1],ib++);
 370   }
 371
 372   TString selectEP = "";
 373   TObjArray resoEParray;
 374   resoEParray.SetOwner();
 375
 376   TString outFilename = "";
 377   TObjArray* optArr = furtherOpt.Tokenize(" ");
 378   TString currName = "";
 379   TString resoFilename = "";
 380   for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
 381     currName = optArr->At(iopt)->GetName();
 382     if ( currName.Contains(".root") ) outFilename = currName;
 383     else if ( currName.Contains(".txt") ) resoFilename = currName;
 384     else {
 385       for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
 386         if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
 387       }
 388     }
 389   }
 390   delete optArr;
 391
 392   if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
 393
 394   furtherOpt.ToUpper();
 395
 396   Bool_t v2only = furtherOpt.Contains("V2ONLY");
 397
 398   TAxis* customBins = 0x0;
 399
 400   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
 401   //Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
 402   //Double_t myBins[] = {6., 8., 10.};
 403   Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
 404   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
 405   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
 406   Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
 407   customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
 408
 409   TCanvas *can = 0x0, *showCan = 0x0;
 410   Int_t xshift = 100;
 411   Int_t yshift = 20;
 412   Int_t igroup1 = 0;
 413   Int_t igroup2 = 0;
 414
 415   TList outList;
 416
 417   //
 418   /// Read plane resolution file (if present)
 419   //
 420   TObjArray resoCentrality;
 421   resoCentrality.SetOwner();
 422   TArrayD resolution[3];
 423   for ( Int_t ires=0; ires<3; ires++ ) {
 424     resolution[ires].Set(50);
 425   }
 426
 427   Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
 428   Bool_t internalReso = kTRUE;
 429   if ( externalReso ) {
 430     // If plane resolution file exists, fill resolution
 431     TString currLine = "";
 432     std::ifstream inFile(resoFilename.Data());
 433     if (inFile.is_open()) {
 434       ib = 0;
 435       while (! inFile.eof() ) {
 436         currLine.ReadLine(inFile,kFALSE);
 437         if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
 438         TObjArray* array = currLine.Tokenize(" ");
 439         resoCentrality.Add(new TObjString(array->At(0)->GetName()));
 440         for ( Int_t ires=0; ires<3; ires++ ) {
 441           resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
 442         }
 443         delete array;
 444         ib++;
 445       } // ! EOF
 446     } // file is open
 447     for ( Int_t ires=0; ires<3; ires++ ) {
 448       resolution[ires].Set(ib);
 449     }
 450   } // file exists
 451   else {
 452     // If no plane resolution file exist,
 453     // use the centralities in terminateOption
 454     // and set the resolution to 1 and errors to 0.
 455
 456     for ( Int_t ires=0; ires<3; ires++ ) {
 457       resolution[ires].Set(centralityClasses->GetEntries());
 458       Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
 459       resolution[ires].Reset(resetVal);
 460     }
 461     TArrayD currCos(3);
 462     for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 463       resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
 464       Int_t icos = 0;
 465       Double_t sumRelErrSquare = 0.;
 466       TString hResoTitle = "";
 467       for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
 468         TString resoDet1 = resoEParray.At(isub)->GetName();
 469         for ( Int_t jsub=isub+1; jsub<resoEParray.GetEntries(); jsub++ ) {
 470           TString resoDet2 = resoEParray.At(jsub)->GetName();
 471           TH1* histo = 0x0;
 472           // Try both combinations: we want the first two values to contain the selectedEP!
 473           TString resoDet = "";
 474           for ( Int_t icomb=0; icomb<2; icomb++ ) {
 475             TString hName = "hResoSub3";
 476             resoDet = ( icomb == 0 ) ? resoDet1 + resoDet2 : resoDet2 + resoDet1;
 477             hName += resoDet;
 478             histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),hName));
 479             if ( histo ) break;
 480           }
 481           if ( ! histo || histo->Integral() == 0 ) continue;
 482           currCos[icos] = histo->GetMean();
 483           if ( currCos[icos] <= 0. ) continue;
 484           Double_t relErr = histo->GetMeanError() / currCos[icos];
 485           sumRelErrSquare += relErr * relErr;
 486           icos++;
 487           hResoTitle += Form("%s ", resoDet.Data());
 488         } // loop on sub events
 489       } // loop on sub events
 490       if ( icos != 3  ) {
 491         printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
 492         internalReso = kFALSE;
 493         continue;
 494       }
 495       resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
 496       resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
 497       printf("Resolution (%s) %s  %g +- %g\n", hResoTitle.Data(), centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
 498     }
 499   }
 500
 501   Bool_t hasResolution = ( externalReso || internalReso );
 502
 503
 504   TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
 505   ib = 0;
 506   for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
 507     TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
 508     if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
 509     resoCentrArray.AddAt(range[1], ib++);
 510   }
 511
 512   const Int_t kNresoCentr = resoCentrality.GetEntries();
 513
 514   //
 515   // Create gridSparse array
 516   //
 517   AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
 518   TString stepName[2];
 519   for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
 520     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 521       gridSparseArray[istep*kNresoCentr+icent] = 0x0;
 522     }
 523
 524     AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
 525     if ( ! cfContainer ) continue;
 526     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 527       stepName[istep] = cfContainer->GetStepTitle(istep);
 528       gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
 529       if ( ! customBins ) customBins = gridSparseArray[istep*kNresoCentr+icent]->GetAxis(kHvarPt);
 530     } // loop on step
 531   } // loop on centrality
 532
 533
 534   Bool_t isMC = furtherOpt.Contains("MC");
 535   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
 536   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
 537   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
 538   if ( ! isMC ) srcColors[kUnidentified] = 1;
 539
 540
 541   //
 542   // Get histograms
 543   //
 544   const Int_t kNcentralities = centralityClasses->GetEntries();
 545   const Int_t kNhistos = kNsteps*kNtrackSources*kNcentralities;
 546   TObjArray ptVsDeltaPhiList(kNhistos);
 547   TObjArray ptVsPhiList(kNhistos);
 548   TArrayD wgtReso[3];
 549   for ( Int_t ires=0; ires<3; ires++ ) {
 550     wgtReso[ires].Set(kNhistos);
 551     wgtReso[ires].Reset(0.);
 552   }
 553   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 554     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 555       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
 556       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 557         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
 558         TH1* hPtVsDeltaPhi = 0x0;
 559         TH1* hPtVsPhi = 0x0;
 560         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
 561         Double_t nMuons = 0.;
 562         TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
 563         for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
 564           if ( resoCentrArray[rcent] < centralityArray[icent] ||
 565               resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
 566           AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
 567           if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
 568
 569           SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
 570           SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
 571           TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
 572
 573           Double_t currYield = auxHisto->Integral();
 574           if ( currYield == 0. ) continue;
 575           nMuons += currYield;
 576           debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
 577           for ( Int_t ires=0; ires<3; ires++ ) {
 578             wgtReso[ires][idx] += resolution[ires][rcent] * currYield;
 579             debugString += Form("  %g", resolution[ires][rcent]);
 580           }
 581           if ( ! hPtVsDeltaPhi ) hPtVsDeltaPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsDeltaPhi_%s", baseNameCent.Data())));
 582           else hPtVsDeltaPhi->Add(auxHisto);
 583           delete auxHisto;
 584
 585           auxHisto = gridSparse->Project(kHvarPhi,kHvarPt);
 586           if ( ! hPtVsPhi ) hPtVsPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsPhi_%s", baseNameCent.Data())));
 587           else hPtVsPhi->Add(auxHisto);
 588           delete auxHisto;
 589
 590         } // loop on resolution centralities
 591
 592         if ( nMuons == 0. ) continue;
 593         debugString += Form("\nWeighted reso ");
 594         for ( Int_t ires=0; ires<3; ires++ ) {
 595           wgtReso[ires][idx] = wgtReso[ires][idx]/nMuons;
 596           debugString += Form("  %g", wgtReso[ires][idx]);
 597         }
 598         AliDebug(1,debugString.Data());
 599
 600         ptVsDeltaPhiList.AddAt(hPtVsDeltaPhi,idx);
 601         ptVsPhiList.AddAt(hPtVsPhi,idx);
 602       } // loop on centralities
 603     } // loop on sources
 604   } // loop on steps
 605
 606
 607
 608   /////////////////////
 609   // Plot resolution //
 610   /////////////////////
 611   if ( hasResolution ) {
 612     currName = Form("reso_%s", selectEP.Data());
 613     currName += externalReso ? Form("_external") : Form("_sub_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
 614     TGraphErrors* resoGraph = new TGraphErrors(kNresoCentr);
 615     resoGraph->SetName(currName.Data());
 616     currName.Append("_can");
 617     can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 618     can->cd();
 619
 620     for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
 621       resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
 622       resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
 623     }
 624     resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
 625     resoGraph->GetXaxis()->SetTitle("Centrality");
 626     resoGraph->GetYaxis()->SetTitle("Resolution");
 627     resoGraph->Draw("apz");
 628     outList.Add(resoGraph);
 629
 630     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 631       for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 632         currName = Form("wgtReso_%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
 633         resoGraph = new TGraphErrors();
 634         resoGraph->SetName(currName.Data());
 635         for ( Int_t icent=0; icent<kNcentralities; icent++ ) {
 636           Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
 637           if ( wgtReso[0][idx] == 0. ) continue;
 638           resoGraph->SetPoint(icent, 0.5*(centralityArray[icent+1] + centralityArray[icent]), wgtReso[0][idx]);
 639           resoGraph->SetPointError(icent, 0.5*(centralityArray[icent+1] - centralityArray[icent]), wgtReso[1][idx]);
 640         } // loop on centralities
 641         if ( resoGraph->GetN() == 0 ) continue;
 642         printf("New reso %s\n", resoGraph->GetName());
 643         resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
 644         resoGraph->GetXaxis()->SetTitle("Centrality");
 645         resoGraph->GetYaxis()->SetTitle("Resolution");
 646         resoGraph->Draw("pz");
 647         outList.Add(resoGraph);
 648       } // loop on sources
 649     } // loop on steps
 650   }
 651
 652
 653
 654   TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
 655   Int_t drawOrder[3] = {0, 2, 1};
 656   TString drawOpt[3] = {"apz","pz","a2"};
 657   Int_t nTypes = ( hasResolution ) ? 3 : 1;
 658
 659   TString histoName = "", histoPattern = "";
 660   ///////////
 661   // Flow  //
 662   ///////////
 663
 664   gStyle->SetOptFit(1111);
 665   TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))";
 666   //
 667   TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
 668   if  ( v2only ) func->SetParNames("scale", "v2");
 669   else func->SetParNames("scale","v1", "v2", "v3");
 670   Int_t v2par = ( v2only ) ? 1 : 2;
 671
 672   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 673     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 674       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
 675       TGraphErrors* flowVsCentrality[3];
 676       for ( Int_t itype=0; itype<nTypes; itype++ ) {
 677         histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
 678         flowVsCentrality[itype] = new TGraphErrors();
 679         flowVsCentrality[itype]->SetName(histoName.Data());
 680         flowVsCentrality[itype]->SetTitle(histoName.Data());
 681       }
 682       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 683         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
 684         TH2* histo2D = static_cast<TH2*> (ptVsDeltaPhiList.At(idx));
 685         if ( ! histo2D ) continue;
 686         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
 687
 688         TAxis* ptAxis = histo2D->GetYaxis();
 689
 690         Int_t ipad = 0;
 691         TGraphErrors* flowVsPt[3];
 692         for ( Int_t itype=0; itype<nTypes; itype++ ) {
 693           histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
 694           flowVsPt[itype] = new TGraphErrors();
 695           flowVsPt[itype]->SetName(histoName.Data());
 696           flowVsPt[itype]->SetTitle(histoName.Data());
 697         }
 698
 699         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
 700           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
 701           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
 702           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
 703           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
 704           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
 705           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
 706           TH1* projHisto = histo2D->ProjectionX(Form("%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
 707           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
 708           if ( projHisto->Integral() < 50 ) break;
 709           if ( ipad%4 == 0 ) {
 710             currName = projHisto->GetName();
 711             currName.Append("_can");
 712             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 713             showCan->Divide(2,2);
 714             ipad = 0;
 715           }
 716           showCan->cd(++ipad);
 717           if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
 718           else {
 719             func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
 720             func->FixParameter(1,0.);
 721           }
 722           projHisto->Fit(func,"R");
 723           for ( Int_t itype=0; itype<nTypes; itype++ ) {
 724             TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
 725             Double_t resoVal = ( itype == 0 ) ? 1. : wgtReso[0][idx];
 726             Double_t resoErr = ( itype == 0 ) ? 0. : wgtReso[itype][idx];
 727             Double_t rawVal = func->GetParameter(v2par);
 728             Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
 729             Double_t yVal = rawVal/resoVal;
 730             Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
 731             Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
 732             Int_t ipoint = currGraph->GetN();
 733             currGraph->SetPoint(ipoint,xVal,yVal);
 734             if ( itype > 0 && ipt != 0 ) {
 735               // For the plot vs. pt the resolution error is fully correlated.
 736               // Hence we do not use add it bin by bin, but rather write it
 737               // on the title
 738               resoErr = 0.;
 739               TString graphTitle = currGraph->GetTitle();
 740               if ( ! graphTitle.Contains("|") ) {
 741                 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", wgtReso[1][idx]/wgtReso[0][idx], wgtReso[2][idx]/wgtReso[0][idx]);
 742                 currGraph->SetTitle(graphTitle.Data());
 743               }
 744             }
 745             Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
 746             Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
 747             Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
 748             currGraph->SetPointError(ipoint,xErr,yErr);
 749           } // loop on type
 750         } // loop on pt bins
 751         for ( Int_t itype=0; itype<nTypes; itype++ ) {
 752           Int_t currType = drawOrder[itype];
 753           if ( itype < 2 ) {
 754             currName = flowVsPt[currType]->GetName();
 755             currName.Append("Can");
 756             can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 757             igroup2++;
 758             can->cd();
 759           }
 760           flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
 761           flowVsPt[currType]->GetYaxis()->SetTitle("v2");
 762           flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
 763           flowVsPt[currType]->SetFillStyle(0);
 764           flowVsPt[currType]->Draw(drawOpt[currType].Data());
 765           outList.Add(flowVsPt[currType]);
 766         } // loop on types
 767       } // loop on centrality
 768       for ( Int_t itype=0; itype<nTypes; itype++ ) {
 769         Int_t currType = drawOrder[itype];
 770         if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
 771         if ( itype < 2 ) {
 772           currName = flowVsCentrality[currType]->GetName();
 773           currName.Append("Can");
 774           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 775           can->cd();
 776         }
 777         flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
 778         flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
 779         flowVsCentrality[currType]->SetFillStyle(0);
 780         flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
 781         flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
 782         outList.Add(flowVsCentrality[currType]);
 783         igroup2++;
 784       } // loop on type
 785       igroup1++;
 786       igroup2 = 0;
 787     } // loop on track sources
 788   } // loop on steps
 789
 790   ///////////////////////
 791   // Event plane check //
 792   ///////////////////////
 793   igroup1++;
 794   igroup2 = 0;
 795   Int_t icolor=0;
 796   currName ="eventPlaneDistrib";
 797   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 798   can->SetGridy();
 799   igroup2++;
 800
 801   const Int_t kNfits = 2;
 802   TString fitTypeName[kNfits] = {"Cos","Sin"};
 803
 804   Bool_t addOffsetToCentrality = kFALSE;
 805   Double_t offsetStep = 0.1;
 806
 807   TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
 808   leg->SetBorderSize(1);
 809
 810   TGraphErrors* epCheck[kNfits];
 811   for ( Int_t itype=0; itype<kNfits; itype++ ) {
 812     histoName = Form("checkFourier_%s",fitTypeName[itype].Data());
 813     epCheck[itype] = new TGraphErrors();
 814     epCheck[itype]->SetName(histoName.Data());
 815 //    epCheck[itype]->SetTitle(histoName.Data());
 816   }
 817   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 818     TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
 819     if ( ! histo ) continue;
 820     histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
 821     if ( histo->Integral() < 50. ) continue;
 822     if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
 823
 824     for ( Int_t itype=0; itype<kNfits; itype++ ) {
 825       TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
 826       TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), histo->GetXaxis()->GetBinLowEdge(1), histo->GetXaxis()->GetBinUpEdge(histo->GetXaxis()->GetNbins()));
 827       checkFunc->SetParameters(histo->Integral(), 0.);
 828       histo->Fit(checkFunc,"R0");
 829       TGraphErrors* currGraph = epCheck[itype];
 830       Double_t yVal = checkFunc->GetParameter(1);
 831       Double_t yErr = checkFunc->GetParError(1);
 832       Double_t xVal = 0.5*(centralityArray[icent+1] + centralityArray[icent]);
 833       Double_t xErr = 0.5*(centralityArray[icent+1] - centralityArray[icent]);
 834       Int_t ipoint = currGraph->GetN();
 835       currGraph->SetPoint(ipoint,xVal,yVal);
 836       currGraph->SetPointError(ipoint,xErr,yErr);
 837     } // loop on type
 838     histo->Fit("pol0","R0");
 839     TF1* pol0 = histo->GetFunction("pol0");
 840     histo->Scale(1./pol0->GetParameter(0));
 841     Double_t offset = ( addOffsetToCentrality ) ? offsetStep*(Double_t)icolor : 0.;
 842     TGraphErrors* graph = new TGraphErrors();
 843     graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
 844     for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
 845       Int_t ibin = ipoint+1;
 846       graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
 847       graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
 848     }
 849     graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
 850     graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
 851     TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
 852     if ( addOffsetToCentrality ) {
 853       graph->GetYaxis()->SetRangeUser(1.-offsetStep, 1.+1.5*offsetStep*(Double_t)centralityClasses->GetEntries());
 854       graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
 855       legTitle += Form(" ( + %g)", offset);
 856     }
 857     Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
 858     graph->SetLineColor(currColor);
 859     graph->SetMarkerColor(currColor);
 860     graph->SetMarkerStyle(20+icolor);
 861     TString currDrawOpt = "pz";
 862     if ( icolor == 0 ) currDrawOpt += "a";
 863     graph->Draw(currDrawOpt.Data());
 864     legTitle.ReplaceAll("_"," - ");
 865     leg->AddEntry(graph,legTitle.Data(),"lp");
 866     if ( addOffsetToCentrality ) {
 867       TLatex* latex = new TLatex(histo->GetXaxis()->GetXmax()+0.5*histo->GetXaxis()->GetBinWidth(1), 1.+(Double_t)offset, Form("#chi^{2}/NDF = %.3g", pol0->GetChisquare()/(Double_t)pol0->GetNDF()));
 868       latex->SetTextColor(currColor);
 869       latex->SetTextSize(0.025);
 870       latex->Draw("same");
 871     }
 872     icolor++;
 873   }
 874   leg->Draw("same");
 875
 876   currName = "epCheckCan";
 877   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 878   can->SetGridy();
 879   igroup2++;
 880   leg = new TLegend(0.7,0.7,0.95,0.95);
 881   leg->SetBorderSize(1);
 882   for ( Int_t itype=0; itype<kNfits; itype++ ) {
 883     epCheck[itype]->SetLineColor(srcColors[itype]);
 884     epCheck[itype]->SetMarkerColor(srcColors[itype]);
 885     epCheck[itype]->SetMarkerStyle(20+itype);
 886     epCheck[itype]->GetXaxis()->SetTitle("Centrality (%)");
 887     TString currDrawOpt = "pz";
 888     if ( can->GetListOfPrimitives()->GetEntries() == 0 ) currDrawOpt.Prepend("a");
 889     epCheck[itype]->Draw(currDrawOpt.Data());
 890     leg->AddEntry(epCheck[itype],Form("<%s(2 #psi_{%s})>", fitTypeName[itype].Data(), selectEP.Data()),"lp");
 891   }
 892   leg->Draw("same");
 893
 894   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 895     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 896       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
 897       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 898         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
 899         TH2* histo2D = static_cast<TH2*> (ptVsPhiList.At(idx));
 900         if ( ! histo2D ) continue;
 901         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
 902
 903         TAxis* ptAxis = histo2D->GetYaxis();
 904
 905         Int_t ipad = 0;
 906         TGraphErrors* detEffect[kNfits];
 907         for ( Int_t itype=0; itype<kNfits; itype++ ) {
 908           histoName = Form("checkResoVsPt_%s_%s", baseNameCent.Data(), fitTypeName[itype].Data());
 909           detEffect[itype] = new TGraphErrors();
 910           detEffect[itype]->SetName(histoName.Data());
 911           detEffect[itype]->SetTitle(histoName.Data());
 912         }
 913
 914         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
 915           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
 916           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
 917           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
 918           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
 919           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
 920           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
 921           TH1* projHisto = histo2D->ProjectionX(Form("checkReso_%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
 922           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
 923           if ( projHisto->Integral() < 50 ) break;
 924           if ( ipad%4 == 0 ) {
 925             currName = projHisto->GetName();
 926             currName.Append("_can");
 927             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 928             showCan->Divide(2,2);
 929             ipad = 0;
 930           }
 931           showCan->cd(++ipad);
 932           for ( Int_t itype=0; itype<kNfits; itype++ ) {
 933             TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
 934             TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), projHisto->GetXaxis()->GetBinLowEdge(1), projHisto->GetXaxis()->GetBinUpEdge(projHisto->GetXaxis()->GetNbins()));
 935             checkFunc->SetParameters(projHisto->Integral(), 0.);
 936             projHisto->Fit(checkFunc,"R");
 937             TGraphErrors* currGraph = detEffect[itype];
 938             Double_t yVal = checkFunc->GetParameter(1);
 939             Double_t yErr = checkFunc->GetParError(1);
 940             Double_t xVal = 0.5*(ptMax + ptMin);
 941             Double_t xErr = 0.5*(ptMax - ptMin);
 942             Int_t ipoint = currGraph->GetN();
 943             currGraph->SetPoint(ipoint,xVal,yVal);
 944             currGraph->SetPointError(ipoint,xErr,yErr);
 945           } // loop on type
 946         } // loop on pt bins
 947         for ( Int_t itype=0; itype<kNfits; itype++ ) {
 948           currName = detEffect[itype]->GetName();
 949           currName.Append("Can");
 950           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 951           igroup2++;
 952           can->cd();
 953           detEffect[itype]->GetXaxis()->SetTitle(ptAxis->GetTitle());
 954           detEffect[itype]->GetYaxis()->SetTitle("v2");
 955 //          detEffect[itype]->GetYaxis()->SetRangeUser(0., 0.25);
 956           detEffect[itype]->SetFillStyle(0);
 957           detEffect[itype]->Draw("ap");
 958 //          outList.Add(detEffect[itype]);
 959         } // loop on types
 960       } // loop on centrality
 961     } // loop on track sources
 962   } // loop on steps
 963
 964   delete centralityClasses;
 965
 966   //////////////////////
 967   // Event statistics //
 968   //////////////////////
 969   printf("\nTotal analyzed events:\n");
 970   TString evtSel = Form("trigger:%s", trigClassName.Data());
 971   fEventCounters->PrintSum(evtSel.Data());
 972   printf("Physics selected analyzed events:\n");
 973   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
 974   fEventCounters->PrintSum(evtSel.Data());
 975
 976   TString countPhysSel = "any";
 977   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
 978   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
 979   countPhysSel.Prepend("selected:");
 980   printf("Analyzed events vs. centrality:\n");
 981   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
 982   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
 983
 984   if ( ! outFilename.IsNull() ) {
 985     printf("\nWriting output file %s\n", outFilename.Data());
 986     TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
 987     outList.Write();
 988     file->Close();
 989   }
 990
 991   delete customBins;
 992 }