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 "AliMuonTrackCuts.h"
  80
  81
  82 /// \cond CLASSIMP
  83 ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
  84 /// \endcond
  85
  86 using std::ifstream;
  87
  88 //________________________________________________________________________
  89 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
  90   AliVAnalysisMuon(),
  91   fEPKeys(0x0),
  92   fRandom(0x0)
  93 {
  94   /// Default ctor.
  95 }
  96
  97 //________________________________________________________________________
  98 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
  99   AliVAnalysisMuon(name, cuts),
 100   fEPKeys(0x0),
 101   fRandom(0x0)
 102 {
 103   //
 104   /// Constructor.
 105   //
 106   TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
 107   fEPKeys = epMethods.Tokenize(" ");
 108 }
 109
 110
 111 //________________________________________________________________________
 112 AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
 113 {
 114   //
 115   /// Destructor
 116   //
 117   delete fEPKeys;
 118   delete fRandom;
 119 }
 120
 121 //___________________________________________________________________________
 122 void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
 123 {
 124   //
 125   /// Create outputs
 126   //
 127
 128   if ( ! fRandom ) fRandom = new TRandom3();
 129   // Set seed here, to avoid having the same seed for all jobs in grid
 130   fRandom->SetSeed(0);
 131
 132   Int_t nPtBins = 160;
 133   Double_t ptMin = 0., ptMax = 80.;
 134   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
 135
 136   Int_t nEtaBins = 25;
 137   Double_t etaMin = -4.5, etaMax = -2.;
 138   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
 139
 140   Int_t nPhiBins = 36;
 141   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
 142   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
 143
 144   Int_t nDeltaPhiBins = 18;
 145   Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
 146   TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
 147
 148   Int_t nChargeBins = 2;
 149   Double_t chargeMin = -2, chargeMax = 2.;
 150   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
 151
 152   Int_t nMotherTypeBins = kNtrackSources;
 153   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
 154   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
 155
 156   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
 157   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
 158   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
 159   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
 160   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
 161
 162   Int_t iobj=0;
 163   TH1* histo = 0x0;
 164   TString currTitle = "";
 165   for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
 166     TString currEP = fEPKeys->At(iep)->GetName();
 167     AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
 168
 169     TString currName = "";
 170     for ( Int_t idim = 0; idim<kNvars; idim++){
 171       currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
 172       currName.ReplaceAll("()","");
 173
 174       cfContainer->SetVarTitle(idim, currName.Data());
 175       cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
 176     }
 177
 178     TString stepTitle[kNsteps] = {"reconstructed", "generated"};
 179
 180     TAxis* currAxis = 0x0;
 181     for (Int_t istep=0; istep<kNsteps; istep++){
 182       cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
 183       AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
 184
 185       currAxis = gridSparse->GetAxis(kHvarMotherType);
 186       for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
 187         currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
 188       }
 189     }
 190
 191     AddObjectToCollection(cfContainer, iobj++);
 192
 193     Bool_t isTPC = currEP.Contains("TPC");
 194     Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
 195     Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
 196
 197     histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
 198     histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
 199     AddObjectToCollection(histo, iobj++);
 200
 201     for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
 202       TString auxEP = fEPKeys->At(jep)->GetName();
 203       currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
 204       histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
 205       histo->GetXaxis()->SetTitle(currTitle.Data());
 206       AddObjectToCollection(histo, iobj++);
 207     }
 208
 209     currTitle = Form("cos(#psi_{1} - #psi_{2})");
 210     histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
 211     histo->GetXaxis()->SetTitle(currTitle.Data());
 212
 213     AddObjectToCollection(histo, iobj++);
 214   } // loop on EPs
 215
 216   fMuonTrackCuts->Print("mask");
 217 }
 218
 219 //________________________________________________________________________
 220 void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
 221 {
 222   //
 223   /// Fill output objects
 224   //
 225
 226   //
 227   // Base event cuts
 228   //
 229   if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
 230   //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
 231   //if ( isPileupFromSPD ) return;
 232   if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) 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     } // loop on trigger
 279   } // loop on EP
 280
 281   //
 282   // Fill track info
 283   //
 284   Double_t containerInput[kNvars];
 285   AliVParticle* track = 0x0;
 286   for ( Int_t istep = 0; istep<2; ++istep ) {
 287     Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
 288     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
 289       track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
 290
 291       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
 292       if ( ! isSelected ) continue;
 293
 294       Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
 295
 296
 297       containerInput[kHvarPt]         = track->Pt();
 298       containerInput[kHvarEta]        = track->Eta();
 299       containerInput[kHvarPhi]        = track->Phi();
 300       containerInput[kHvarCharge]     = track->Charge()/3.;
 301       containerInput[kHvarMotherType] = (Double_t)trackSrc;
 302
 303       for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
 304         if ( psiPlane[iep] < kLowestVal ) continue;
 305         TString currEP = fEPKeys->At(iep)->GetName();
 306         Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
 307
 308         // The difference between phi - psiPlane must be between (0, pi)
 309         // Since the reaction plain is symmetric in pi,
 310         // we can do in such a way that:
 311         // -pi < delta_phi + N*pi < pi
 312         // We choose N accrodingly
 313         Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
 314         deltaPhi += nPi * TMath::Pi();
 315         containerInput[kHvarDeltaPhi] = deltaPhi;
 316
 317         for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 318           TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 319           if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
 320           ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
 321         } // loop on selected trigger classes
 322       } // loop on event plane
 323     } // loop on tracks
 324   } // loop on container steps
 325 }
 326
 327
 328 //________________________________________________________________________
 329 TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
 330 {
 331   //
 332   /// Get the range from string
 333   //
 334   TArrayD centralityRange(2);
 335   centralityRange.Reset(-999.);
 336   TObjArray* array = sRange.Tokenize("_");
 337   if ( array->GetEntries() == 2 ) {
 338     for ( Int_t ientry=0; ientry<2; ientry++ ) {
 339       centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
 340     }
 341   }
 342   delete array;
 343   return centralityRange;
 344 }
 345
 346 //________________________________________________________________________
 347 void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
 348   //
 349   /// Draw some histograms at the end.
 350   //
 351
 352   AliVAnalysisMuon::Terminate("");
 353
 354   if ( ! fMergeableCollection ) return;
 355
 356   TString physSel = fTerminateOptions->At(0)->GetName();
 357   TString trigClassName = fTerminateOptions->At(1)->GetName();
 358   TString centralityRange = fTerminateOptions->At(2)->GetName();
 359   TString furtherOpt = fTerminateOptions->At(3)->GetName();
 360
 361   TObjArray* centralityClasses = centralityRange.Tokenize(",");
 362   TArrayD centralityArray(centralityClasses->GetEntries()+1);
 363   Int_t ib = 0;
 364   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 365     TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
 366     if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
 367     centralityArray.AddAt(range[1],ib++);
 368   }
 369
 370   TString selectEP = "";
 371   TObjArray resoEParray;
 372   resoEParray.SetOwner();
 373
 374   TString outFilename = "";
 375   TObjArray* optArr = furtherOpt.Tokenize(" ");
 376   TString currName = "";
 377   TString resoFilename = "";
 378   for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
 379     currName = optArr->At(iopt)->GetName();
 380     if ( currName.Contains(".root") ) outFilename = currName;
 381     else if ( currName.Contains(".txt") ) resoFilename = currName;
 382     else {
 383       for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
 384         if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
 385       }
 386     }
 387   }
 388   delete optArr;
 389
 390   if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
 391
 392   furtherOpt.ToUpper();
 393
 394   Bool_t v2only = furtherOpt.Contains("V2ONLY");
 395
 396   TAxis* customBins = 0x0;
 397
 398   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
 399   Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
 400   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
 401   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
 402   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
 403   Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
 404   customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
 405
 406   TCanvas *can = 0x0, *showCan = 0x0;
 407   Int_t xshift = 100;
 408   Int_t yshift = 20;
 409   Int_t igroup1 = 0;
 410   Int_t igroup2 = 0;
 411
 412   //
 413   /// Read plane resolution file (if present)
 414   //
 415   TObjArray resoCentrality;
 416   resoCentrality.SetOwner();
 417   TArrayD resolution[3];
 418   for ( Int_t ires=0; ires<3; ires++ ) {
 419     resolution[ires].Set(50);
 420   }
 421
 422   Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
 423   Bool_t internalReso = kTRUE;
 424   if ( externalReso ) {
 425     // If plane resolution file exists, fill resolution
 426     TString currLine = "";
 427     std::ifstream inFile(resoFilename.Data());
 428     if (inFile.is_open()) {
 429       ib = 0;
 430       while (! inFile.eof() ) {
 431         currLine.ReadLine(inFile,kFALSE);
 432         if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
 433         TObjArray* array = currLine.Tokenize(" ");
 434         resoCentrality.Add(new TObjString(array->At(0)->GetName()));
 435         for ( Int_t ires=0; ires<3; ires++ ) {
 436           resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
 437         }
 438         delete array;
 439         ib++;
 440       } // ! EOF
 441     } // file is open
 442     for ( Int_t ires=0; ires<3; ires++ ) {
 443       resolution[ires].Set(ib);
 444     }
 445   } // file exists
 446   else {
 447     // If no plane resolution file exist,
 448     // use the centralities in terminateOption
 449     // and set the resolution to 1 and errors to 0.
 450
 451     for ( Int_t ires=0; ires<3; ires++ ) {
 452       resolution[ires].Set(centralityClasses->GetEntries());
 453       Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
 454       resolution[ires].Reset(resetVal);
 455     }
 456     TArrayD currCos(3);
 457     for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 458       resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
 459       Int_t icos = 0;
 460       Double_t sumRelErrSquare = 0.;
 461       for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
 462         for ( Int_t jsub=0; jsub<resoEParray.GetEntries(); jsub++ ) {
 463           TH1* histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName())));
 464           if ( ! histo || histo->Integral() == 0 ) continue;
 465           currCos[icos] = histo->GetMean();
 466           if ( currCos[icos] <= 0. ) continue;
 467           Double_t relErr = histo->GetMeanError() / currCos[icos];
 468           sumRelErrSquare += relErr * relErr;
 469           icos++;
 470         } // loop on sub events
 471       } // loop on sub events
 472       if ( icos != 3  ) {
 473         printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
 474         internalReso = kFALSE;
 475         continue;
 476       }
 477       resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
 478       resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
 479       printf("Resolution %s  %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
 480     }
 481   }
 482
 483   Bool_t hasResolution = ( externalReso || internalReso );
 484
 485
 486   TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
 487   ib = 0;
 488   for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
 489     TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
 490     if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
 491     resoCentrArray.AddAt(range[1], ib++);
 492   }
 493
 494   if ( hasResolution ) {
 495     currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
 496     can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 497     TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries());
 498     for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
 499
 500       resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
 501       resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
 502     }
 503     resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
 504     resoGraph->GetXaxis()->SetTitle("Centrality");
 505     resoGraph->GetYaxis()->SetTitle("Resolution");
 506     resoGraph->Draw("apz");
 507   }
 508
 509   const Int_t kNresoCentr = resoCentrality.GetEntries();
 510
 511   //
 512   // Create gridSparse array
 513   //
 514   AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
 515   TString stepName[2];
 516   for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
 517     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 518       gridSparseArray[istep*kNresoCentr+icent] = 0x0;
 519     }
 520
 521     AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
 522     if ( ! cfContainer ) continue;
 523     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 524       stepName[istep] = cfContainer->GetStepTitle(istep);
 525       gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
 526     } // loop on step
 527   } // loop on centrality
 528
 529   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
 530
 531   Bool_t isMC = furtherOpt.Contains("MC");
 532   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
 533   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
 534   if ( ! isMC ) srcColors[kUnidentified] = 1;
 535
 536   TList outList;
 537
 538   TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
 539   Int_t drawOrder[3] = {0, 2, 1};
 540   TString drawOpt[3] = {"apz","pz","a2"};
 541   Int_t nTypes = ( hasResolution ) ? 3 : 1;
 542
 543   TString histoName = "", histoPattern = "";
 544   ///////////
 545   // Flow  //
 546   ///////////
 547   gStyle->SetOptFit(1111);
 548   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)))";
 549   //
 550   TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
 551   if  ( v2only ) func->SetParNames("scale", "v2");
 552   else func->SetParNames("scale","v1", "v2", "v3");
 553   Int_t v2par = ( v2only ) ? 1 : 2;
 554
 555   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
 556     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 557       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
 558       TGraphErrors* flowVsCentrality[3];
 559       for ( Int_t itype=0; itype<nTypes; itype++ ) {
 560         histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
 561         flowVsCentrality[itype] = new TGraphErrors();
 562         flowVsCentrality[itype]->SetName(histoName.Data());
 563         flowVsCentrality[itype]->SetTitle(histoName.Data());
 564       }
 565       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 566         TH2* hDeltaPhi2D = 0x0;
 567         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
 568         TArrayD weightedReso(3);
 569         weightedReso.Reset(0.);
 570         Double_t nMuons = 0.;
 571         TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
 572         for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
 573           if ( resoCentrArray[rcent] < centralityArray[icent] ||
 574                resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
 575           AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
 576           if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
 577
 578           SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
 579           SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
 580           TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
 581
 582           Double_t currYield = auxHisto->Integral();
 583           if ( currYield == 0. ) continue;
 584           nMuons += currYield;
 585           debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
 586           for ( Int_t ires=0; ires<3; ires++ ) {
 587             weightedReso[ires] += resolution[ires][rcent] * currYield;
 588             debugString += Form("  %g", resolution[ires][rcent]);
 589           }
 590           histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data());
 591           if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast<TH2*>(auxHisto->Clone(histoName.Data()));
 592           else hDeltaPhi2D->Add(auxHisto);
 593           delete auxHisto;
 594         }
 595
 596         if ( ! hDeltaPhi2D ) continue;
 597
 598         debugString += Form("\nWeighted reso ");
 599         for ( Int_t ires=0; ires<3; ires++ ) {
 600           weightedReso[ires] = weightedReso[ires]/nMuons;
 601           debugString += Form("  %g", weightedReso[ires]);
 602         }
 603         AliDebug(1,debugString.Data());
 604
 605         TAxis* ptAxis = hDeltaPhi2D->GetYaxis();
 606
 607         Int_t ipad = 0;
 608         TGraphErrors* flowVsPt[3];
 609         for ( Int_t itype=0; itype<nTypes; itype++ ) {
 610           histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
 611           flowVsPt[itype] = new TGraphErrors();
 612           flowVsPt[itype]->SetName(histoName.Data());
 613           flowVsPt[itype]->SetTitle(histoName.Data());
 614         }
 615
 616         if ( ! customBins ) customBins = static_cast<TAxis*>(ptAxis->Clone());
 617
 618
 619         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
 620           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
 621           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
 622           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
 623           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
 624           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
 625           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
 626           TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin);
 627           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
 628           if ( projHisto->Integral() < 50 ) break;
 629           if ( ipad%4 == 0 ) {
 630             currName = projHisto->GetName();
 631             currName.Append("_can");
 632             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 633             showCan->Divide(2,2);
 634             ipad = 0;
 635           }
 636           showCan->cd(++ipad);
 637           if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
 638           else {
 639             func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
 640             func->FixParameter(1,0.);
 641           }
 642           projHisto->Fit(func,"R");
 643           for ( Int_t itype=0; itype<nTypes; itype++ ) {
 644             TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
 645             Double_t resoVal = ( itype == 0 ) ? 1. : weightedReso[0];
 646             Double_t resoErr = ( itype == 0 ) ? 0. : weightedReso[itype];
 647             Double_t rawVal = func->GetParameter(v2par);
 648             Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
 649             Double_t yVal = rawVal/resoVal;
 650             Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
 651             Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
 652             Int_t ipoint = currGraph->GetN();
 653             currGraph->SetPoint(ipoint,xVal,yVal);
 654             if ( itype > 0 && ipt != 0 ) {
 655               // For the plot vs. pt the resolution error is fully correlated.
 656               // Hence we do not use add it bin by bin, but rather write it
 657               // on the title
 658               resoErr = 0.;
 659               TString graphTitle = currGraph->GetTitle();
 660               if ( ! graphTitle.Contains("|") ) {
 661                 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]);
 662                 currGraph->SetTitle(graphTitle.Data());
 663               }
 664             }
 665             Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
 666             Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
 667             Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
 668             currGraph->SetPointError(ipoint,xErr,yErr);
 669           } // loop on type
 670         } // loop on pt bins
 671         for ( Int_t itype=0; itype<nTypes; itype++ ) {
 672           Int_t currType = drawOrder[itype];
 673           if ( itype < 2 ) {
 674             currName = flowVsPt[currType]->GetName();
 675             currName.Append("Can");
 676             can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 677             igroup2++;
 678             can->cd();
 679           }
 680           flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
 681           flowVsPt[currType]->GetYaxis()->SetTitle("v2");
 682           flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
 683           flowVsPt[currType]->SetFillStyle(0);
 684           flowVsPt[currType]->Draw(drawOpt[currType].Data());
 685           outList.Add(flowVsPt[currType]);
 686         } // loop on types
 687       } // loop on centrality
 688       for ( Int_t itype=0; itype<nTypes; itype++ ) {
 689         Int_t currType = drawOrder[itype];
 690         if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
 691         if ( itype < 2 ) {
 692           currName = flowVsCentrality[currType]->GetName();
 693           currName.Append("Can");
 694           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 695           can->cd();
 696         }
 697         flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
 698         flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
 699         flowVsCentrality[currType]->SetFillStyle(0);
 700         flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
 701         flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
 702         outList.Add(flowVsCentrality[currType]);
 703         igroup2++;
 704       } // loop on type
 705       igroup1++;
 706       igroup2 = 0;
 707     } // loop on track sources
 708   } // loop on steps
 709
 710   delete customBins;
 711
 712
 713   ///////////////////////
 714   // Event plane check //
 715   ///////////////////////
 716   igroup1++;
 717   igroup2 = 0;
 718   Int_t icolor=0;
 719   currName ="eventPlaneDistrib";
 720   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 721   can->SetGridy();
 722   igroup2++;
 723
 724   Bool_t addOffsetToCentrality = kFALSE;
 725
 726   TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
 727   leg->SetBorderSize(1);
 728   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
 729     TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
 730     if ( ! histo ) continue;
 731     histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
 732     if ( histo->Integral() < 50. ) continue;
 733     if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
 734     histo->Fit("pol0","R0");
 735     TF1* pol0 = histo->GetFunction("pol0");
 736     histo->Scale(1./pol0->GetParameter(0));
 737     Int_t offset = ( addOffsetToCentrality ) ? icolor : 0;
 738     TGraphErrors* graph = new TGraphErrors();
 739     graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
 740     for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
 741       Int_t ibin = ipoint+1;
 742       graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
 743       graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
 744     }
 745     graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
 746     graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
 747     TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
 748     if ( addOffsetToCentrality ) {
 749       graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries());
 750       graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
 751       legTitle += Form(" ( + %i)", offset);
 752     }
 753     Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
 754     graph->SetLineColor(currColor);
 755     graph->SetMarkerColor(currColor);
 756     graph->SetMarkerStyle(20+icolor);
 757     TString currDrawOpt = "pz";
 758     if ( icolor == 0 ) currDrawOpt += "a";
 759     graph->Draw(currDrawOpt.Data());
 760     legTitle.ReplaceAll("_"," - ");
 761     leg->AddEntry(graph,legTitle.Data(),"lp");
 762     if ( addOffsetToCentrality ) {
 763       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()));
 764       latex->SetTextColor(currColor);
 765       latex->SetTextSize(0.025);
 766       latex->Draw("same");
 767     }
 768     icolor++;
 769   }
 770   leg->Draw("same");
 771
 772
 773   //////////////////////
 774   // Event statistics //
 775   //////////////////////
 776   printf("\nTotal analyzed events:\n");
 777   TString evtSel = Form("trigger:%s", trigClassName.Data());
 778   fEventCounters->PrintSum(evtSel.Data());
 779   printf("Physics selected analyzed events:\n");
 780   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
 781   fEventCounters->PrintSum(evtSel.Data());
 782
 783   TString countPhysSel = "any";
 784   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
 785   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
 786   countPhysSel.Prepend("selected:");
 787   printf("Analyzed events vs. centrality:\n");
 788   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
 789   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
 790
 791   if ( ! outFilename.IsNull() ) {
 792     printf("\nWriting output file %s\n", outFilename.Data());
 793     TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
 794     outList.Write();
 795     file->Close();
 796   }
 797 }