PWG/muon/AliAnalysisMuMuNch.cxx

   1 #include "AliAnalysisMuMuNch.h"
   2
   3 /**
   4  *
   5  * \ingroup pwg-muon-mumu
   6  *
   7  * \class AliAnalysisMuMuNch
   8  *
   9  * SPD tracklet to Nch analysis
  10  *
  11  * The idea is that this sub-analysis is used first (within the AliAnalysisTaskMuMu sub-framework)
  12  * in order to compute the number of charged particles within an event, so that information
  13  * can be used in subsequent sub-analysis, like the invariant mass or mean pt ones.
  14  *
  15  */
  16
  17 #include "AliInputEventHandler.h"
  18 #include "AliMultiplicity.h"
  19 #include "AliAODTracklets.h"
  20 #include "AliAnalysisMuonUtility.h"
  21 #include "AliLog.h"
  22 #include "AliAODEvent.h"
  23 #include "AliESDEvent.h"
  24 #include "AliESDUtils.h"
  25 #include "TMath.h"
  26 #include "AliAnalysisMuMuCutRegistry.h"
  27 #include "AliAnalysisMuMuEventCutter.h"
  28 #include "AliAnalysisMuMuCutElement.h"
  29 #include "AliAnalysisMuMuBinning.h"
  30 #include "Riostream.h"
  31 #include "TParameter.h"
  32 #include <set>
  33 #include <utility>
  34 #include "AliMergeableCollection.h"
  35 #include "AliMCEvent.h"
  36 #include "AliAODMCParticle.h"
  37 #include "TAxis.h"
  38 #include "TCanvas.h"
  39 #include "TH2F.h"
  40 #include "TF1.h"
  41 #include "TProfile.h"
  42 #include "TObjString.h"
  43 #include "TRandom3.h"
  44
  45 #include "AliAODMCHeader.h"
  46 #include "AliGenEventHeader.h"
  47 #include "AliGenHijingEventHeader.h"
  48 #include "AliGenDPMjetEventHeader.h"
  49 #include "AliGenPythiaEventHeader.h"
  50 #include "AliGenCocktailEventHeader.h"
  51
  52
  53 namespace {
  54
  55   Double_t SPDgeomR(Double_t* x,Double_t* par) // Eta position of the SPD right edge eta as "seen" from a z vertex position
  56   {
  57     // par[0] = radius of SPD layer
  58
  59     Double_t d = x[0];
  60     Double_t z(0.);
  61     Double_t theta(0.);
  62
  63     if( d < 14.09999 )
  64     {
  65       z = 14.1 - d;
  66       theta = TMath::ATan(par[0]/z);
  67     }
  68
  69     else if( d > 14.09999 )
  70     {
  71       z = d - 14.1;
  72       theta = TMath::Pi() - TMath::ATan(par[0]/z);
  73     }
  74
  75     return -TMath::Log(TMath::Tan(theta/2.));
  76   }
  77
  78   Double_t SPDgeomL(Double_t* x,Double_t* par) // Eta position of the SPD left edge eta as "seen" from a z vertex position
  79   {
  80     // par[0] = radius of SPD layer
  81
  82     Double_t d = x[0];
  83     Double_t z(0.);
  84     Double_t theta(0.);
  85
  86     if( d > -14.09999 )
  87     {
  88       z = 14.1 + d;
  89       theta = TMath::Pi() - TMath::ATan(par[0]/z);
  90     }
  91
  92     if( d < -14.09999 )
  93     {
  94       z = -14.1 - d;
  95       theta = TMath::ATan(par[0]/z);
  96     }
  97
  98     return -TMath::Log(TMath::Tan(theta/2.));
  99   }
 100
 101 }
 102
 103 ////_____________________________________________________________________________
 104 //AliAnalysisMuMuNch::AliAnalysisMuMuNch(TH2* spdCorrection, Double_t etaMin, Double_t etaMax
 105 //                                      , Double_t zMin, Double_t zMax,Bool_t disableHistos, Bool_t computeResolution)
 106 //: AliAnalysisMuMuBase(),
 107 //fSPDCorrection(0x0),
 108 //fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
 109 //fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
 110 //fCurrentEvent(0x0),
 111 //fEtaMin(etaMin),
 112 //fEtaMax(etaMax),
 113 //fZMin(zMin),
 114 //fZMax(zMax),
 115 //fResolution(computeResolution)
 116 //{
 117 //  if ( spdCorrection )
 118 //  {
 119 //    fSPDCorrection = static_cast<TH2F*>(spdCorrection->Clone());
 120 //    fSPDCorrection->SetDirectory(0);
 121 //  }
 122 //  DefineSPDAcceptance();
 123 //
 124 //  if ( disableHistos )
 125 //  {
 126 //    DisableHistograms("*");
 127 //  }
 128 //}
 129
 130 //FIXME: First and second constructor may be ambiguous if we do not set all the arguments
 131
 132 //_____________________________________________________________________________
 133 AliAnalysisMuMuNch::AliAnalysisMuMuNch(TH2F* spdCorrection, TProfile* spdMeanCorrection, Double_t meanTrRef, Double_t etaMin, Double_t etaMax
 134                                        , Double_t zMin, Double_t zMax, Bool_t disableHistos, Bool_t computeResolution)
 135 : AliAnalysisMuMuBase(),
 136 fSPDOneOverAccxEff(0x0),
 137 fSPDCorrectionMap(0x0),
 138 fSPDCorrectionList(0x0),
 139 fSPDMeanTracklets(0x0),
 140 fSPDMeanTrackletsCorrToCompare(0x0),
 141 fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
 142 fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
 143 fCurrentEvent(0x0),
 144 fMeanTrRef(meanTrRef),
 145 fEtaMin(etaMin),
 146 fEtaMax(etaMax),
 147 fEtaMinToCompare(0.),
 148 fEtaMaxToCompare(0.),
 149 fetaRange(),
 150 fZMin(zMin),
 151 fZMax(zMax),
 152 fResolution(computeResolution),
 153 frand(0x0),
 154 fGeneratorHeaderClass(new TString("AliGenDPMjetEventHeader"))
 155 {
 156   //FIXME: Add a protection to avoid an etamin or etamax non multiple of the eta bin size
 157
 158   //FIXME: Add fluctuations to the SPD AccxEff correction method
 159
 160   if ( spdCorrection && spdMeanCorrection )
 161   {
 162     AliWarning("Two methods used to correct tracklets: Data-driven and SPDAccxEff corrected tracklets will be available in the event list, but the histograms in this class will be filled just with the AccxEff corrected values");
 163 //    AliFatal("Cannot use 2 different methods for tracklets correction");
 164   }
 165 //  else if ( spdCorrection )
 166   if ( spdCorrection )
 167   {
 168     AliWarning("SPD AccxEff correction method not ready to be used, fluctuations are missing. Do not trust the results");
 169
 170     if ( fMeanTrRef > 0. && !spdMeanCorrection ) AliWarning("Reference mean nof tracklets argument will not be used: SPD AccxEff correction method in use");
 171
 172     fSPDOneOverAccxEff = static_cast<TH2F*>(spdCorrection->Clone());
 173     fSPDOneOverAccxEff->SetDirectory(0);
 174   }
 175 //  else if ( spdMeanCorrection )
 176   if ( spdMeanCorrection )
 177   {
 178     if ( fMeanTrRef < 0. ) AliWarning("Reference mean nof tracklets argument set to -1: Maximum of the spdMeanCorrection will be used");
 179     else AliWarning(Form("Using %f as reference mean nof tracklets for correction",fMeanTrRef));
 180
 181     frand = new TRandom3();
 182     fSPDMeanTracklets = static_cast<TProfile*>(spdMeanCorrection->Clone());
 183     fSPDMeanTracklets->SetDirectory(0);
 184   }
 185
 186   DefineSPDAcceptance();
 187
 188   if ( disableHistos ) // FIXME: Is this really useful? it breaks when setting to ktrue due to non existence of histos in SetEvent(). Answer: It is useful, it will speed up the task when we want to execute only SetEvent and not fill all the multiplicity histos (i.e. doing J/psi vs multiplicity analysis). The problem is that disabling the histos the method CreateHisto() does not create histosm thats why the task breaks in SetEvent, so fix this
 189   {
 190     DisableHistograms("*");
 191   }
 192 }
 193
 194 //_____________________________________________________________________________
 195 AliAnalysisMuMuNch::AliAnalysisMuMuNch(TProfile* spdMeanCorrection, TProfile* spdMeanCorrectionToCompare, Double_t meanTrRef, Double_t etaMin,
 196                                        Double_t etaMax, Double_t zMin, Double_t zMax, Double_t etaMinToCompare, Double_t etaMaxToCompare,Bool_t disableHistos, Bool_t computeResolution)
 197 : AliAnalysisMuMuBase(),
 198 fSPDOneOverAccxEff(0x0),
 199 fSPDCorrectionMap(0x0),
 200 fSPDCorrectionList(0x0),
 201 fSPDMeanTracklets(0x0),
 202 fSPDMeanTrackletsCorrToCompare(0x0),
 203 fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
 204 fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
 205 fCurrentEvent(0x0),
 206 fMeanTrRef(meanTrRef),
 207 fEtaMin(etaMin),
 208 fEtaMax(etaMax),
 209 fEtaMinToCompare(etaMinToCompare),
 210 fEtaMaxToCompare(etaMaxToCompare),
 211 fetaRange(),
 212 fZMin(zMin),
 213 fZMax(zMax),
 214 fResolution(computeResolution),
 215 frand(0x0),
 216 fGeneratorHeaderClass(new TString("AliGenDPMjetEventHeader"))
 217 {
 218   //FIXME: Add a protection to avoid an etamin or etamax non multiple of the eta bin size
 219
 220   /// This construction is designed to compute everything in etaMin < eta < etaMax and correct with spdMeanCorrection (main eta tange and correction), but when filling the histos in FillHistosForEvent is able to compute the number of tracklets in etaMinToCompare < eta < etaMaxToCompare and correct them with spdMeanCorrectionToCompare (secondary eta range and correction) in order to compare N_{tr}^{|eta|< etaMax} vs N_{tr}^{|eta|< etaMaxToCompare}
 221
 222   if ( !spdMeanCorrection || !spdMeanCorrectionToCompare )
 223   {
 224     AliFatal("Need the two corrections to compare. Maybe you are using the wrong constructor");
 225   }
 226   if ( fEtaMinToCompare < fEtaMin || fEtaMaxToCompare > fEtaMax )
 227   {
 228     AliFatal("Cannot select a eta range to compare wider than the main eta range");
 229   }
 230
 231   if ( fMeanTrRef < 0. ) AliWarning("Reference mean nof tracklets argument set to -1: Maximum of the spdMeanCorrection will be used");
 232   else AliWarning(Form("Using %f as reference mean nof tracklets for correction",fMeanTrRef));
 233
 234   frand = new TRandom3();
 235   fSPDMeanTracklets = static_cast<TProfile*>(spdMeanCorrection->Clone());
 236   fSPDMeanTracklets->SetDirectory(0);
 237
 238   fSPDMeanTrackletsCorrToCompare = static_cast<TProfile*>(spdMeanCorrectionToCompare->Clone());
 239   fSPDMeanTrackletsCorrToCompare->SetDirectory(0);
 240
 241   DefineSPDAcceptance();
 242
 243   if ( disableHistos ) // FIXME: Is this really useful? it breaks when setting to ktrue due to non existence of histos in SetEvent()
 244   {
 245     DisableHistograms("*");
 246   }
 247 }
 248
 249 //_____________________________________________________________________________
 250 AliAnalysisMuMuNch::AliAnalysisMuMuNch(TObjArray* spdCorrectionList, Double_t meanTrRef, Double_t etaMin, Double_t etaMax
 251                                        , Double_t zMin, Double_t zMax,Bool_t disableHistos, Bool_t computeResolution)
 252 : AliAnalysisMuMuBase(),
 253 fSPDOneOverAccxEff(0x0),
 254 fSPDCorrectionMap(0x0),
 255 fSPDCorrectionList(0x0),
 256 fSPDMeanTracklets(0x0),
 257 fSPDMeanTrackletsCorrToCompare(0x0),
 258 fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
 259 fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
 260 fCurrentEvent(0x0),
 261 fMeanTrRef(meanTrRef),
 262 fEtaMin(etaMin),
 263 fEtaMax(etaMax),
 264 fEtaMinToCompare(0.),
 265 fEtaMaxToCompare(0.),
 266 fetaRange(),
 267 fZMin(zMin),
 268 fZMax(zMax),
 269 fResolution(computeResolution),
 270 frand(0x0),
 271 fGeneratorHeaderClass(new TString("AliGenDPMjetEventHeader"))
 272 {
 273   //FIXME: Add a protection to avoid an etamin or etamax non multiple of the eta bin size
 274
 275   // Uses a different correction for each group of runs (both SPD AccxEff OR mean tracklets are supported)
 276
 277   if ( spdCorrectionList ) DefineSPDCorrectionMap(spdCorrectionList);
 278   else AliFatal("No SPD correction list provided");
 279
 280   if ( fMeanTrRef < 0. ) AliWarning("Reference mean nof tracklets argument set to -1: Maximum of the spdMeanCorrection will be used");
 281   else AliWarning(Form("Using %f as reference mean nof tracklets for correction",fMeanTrRef));
 282
 283   frand = new TRandom3();
 284
 285   DefineSPDAcceptance();
 286
 287   if ( disableHistos ) // FIXME: Is this really useful? it breaks when setting to ktrue due to non existence of histos in SetEvent()
 288   {
 289     DisableHistograms("*");
 290   }
 291 }
 292
 293 //_____________________________________________________________________________
 294 AliAnalysisMuMuNch::~AliAnalysisMuMuNch()
 295 {
 296   delete fSPDOneOverAccxEff;
 297   delete fSPDCorrectionMap;
 298   delete fSPDCorrectionList;
 299   delete fSPDMeanTracklets;
 300   delete fSPDMeanTrackletsCorrToCompare;
 301   delete fEtaAxis;
 302   delete fZAxis;
 303   delete fSPD1LR;
 304   delete fSPD1LL;
 305   delete fSPD2LR;
 306   delete fSPD2LL;
 307   delete frand;
 308   delete fGeneratorHeaderClass;
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliAnalysisMuMuNch::DefineGeneratorName(const char* genName)
 313 {
 314 //  TString sgenName(genName);
 315 //  if ( sgenNam.Contains("pythia") ) fGeneratorHeaderClass = "AliGenPythiaEventHeader";
 316 //  else if ( sgenNam.Contains("dpmjet") ) fGeneratorHeaderClass = "AliGenDPMjetEventHeader";
 317 //  else if ( sgenNam.Contains("dpmjet") ) fGeneratorHeaderClass = "AliGenHijingEventHeader";
 318
 319   fGeneratorHeaderClass->Form("%s",genName);
 320
 321   std::cout << " Will use " << fGeneratorHeaderClass->Data() << " tracks" << std::endl;
 322 }
 323
 324 //_____________________________________________________________________________
 325 void AliAnalysisMuMuNch::DefineHistogramCollection(const char* eventSelection,
 326                                                    const char* triggerClassName,
 327                                                    const char* centrality)
 328 {
 329   // Define multiplicity histos
 330
 331   if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuNch") )
 332   {
 333     return;
 334   }
 335
 336   // dummy histogram to signal that we already defined all our histograms (see above)
 337   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"AliAnalysisMuMuNch","Dummy semaphore",1,0,1);
 338
 339   Double_t multMin = -0.5;  //Tracklets multiplicity range
 340   Double_t multMax = 500.5;
 341   Int_t nbinsMult = GetNbins(multMin,multMax,1.);
 342
 343   Double_t phimin = 0.; //Phi range
 344   Double_t phimax = 2*TMath::Pi();
 345   Int_t nphibins = GetNbins(phimin,phimax,0.05);
 346
 347   if ( !fSPDMeanTracklets && !fSPDOneOverAccxEff && fResolution ) // Resolution histos
 348   {
 349     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"EtaRes","#eta resolution;#eta_{Reco} - #eta_{MC};Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 350
 351     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiRes","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*2,-phimax,phimax);
 352
 353     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiResShifted","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*4,-phimax/4,phimax/4);
 354
 355     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"EtaResVsZ","#eta resolution vs MC Zvertex;ZVertex (cm);#eta_{Reco} - #eta_{MC}",(fZAxis->GetNbins())*20,fZAxis->GetXmin(),fZAxis->GetXmax(),(fEtaAxis->GetNbins())*8,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 356
 357     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiResVsZ","#phi resolution vs MC Zvertex;ZVertex (cm);#phi_{Reco} - #phi_{MC}",(fZAxis->GetNbins())*20,fZAxis->GetXmin(),fZAxis->GetXmax(),nphibins*4,-phimax/4,phimax/4);
 358
 359     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"EtaResVsnC","#eta resolution vs Nof contributors to SPD vertex;NofContributors;#eta_{Reco} - #eta_{MC}",200,0,200,(fEtaAxis->GetNbins())*8,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 360
 361     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiResVsnC","#phi resolution vs Nof contributors to SPD vertex;NofContributors;#phi_{Reco} - #phi_{MC}",200,0,200,nphibins*4,-phimax/4,phimax/4);
 362
 363     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"SPDZvResVsnC","SPD Zvertex resolution vs Nof contributors;NofContributors;Zvertex_{Reco} - Zvertex_{MC} (cm)",200,0,200,(fZAxis->GetNbins())*20,fZAxis->GetXmin(),fZAxis->GetXmax());
 364
 365     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"SPDZvResVsMCz","SPD Zvertex resolution vs MC z vertex;MC Zvertex (cm);Zvertex_{Reco} - Zvertex_{MC} (cm)",(fZAxis->GetNbins())*20,fZAxis->GetXmin(),fZAxis->GetXmax(),(fZAxis->GetNbins())*10,fZAxis->GetXmin(),fZAxis->GetXmax());
 366
 367     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Phi","Reco #phi distribution;#phi;Counts",nphibins,phimin,phimax);
 368
 369     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCPhi","MC #phi distribution;#phi;Counts",nphibins,phimin,phimax);
 370
 371     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","Reco #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 372
 373     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCEta","MC #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 374
 375     return; // When computing resolutions we don't want to create the rest of histos
 376   }
 377
 378   TString TrackletsCorrectedName("");
 379   TString TrackletsCorrectedAxisName("");
 380   if ( fSPDMeanTracklets && !fSPDOneOverAccxEff )
 381   {
 382     TrackletsCorrectedName = "Number of corrected tracklets";
 383     TrackletsCorrectedAxisName = "N_{tr}^{corr}";
 384   }
 385   else if( fSPDOneOverAccxEff )
 386   {
 387     TrackletsCorrectedName = "Number of charged particles";
 388     TrackletsCorrectedAxisName = "N_{ch}";
 389   }
 390
 391   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi","Number of tracklets vs Z vertex vs #phi;Z vertex;#phi",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),nphibins,phimin,phimax);
 392
 393   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta","Number of tracklets vs ZVertex vs #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 394   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta");
 395
 396   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Tracklets",Form("Number of tracklets in |#eta| < %1.1f distribution;N_{Tracklets};N_{events}",fEtaMax),nbinsMult,multMin,multMax);
 397
 398   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NBkgTrackletsVsZVertexVsEta","Number of background tracklets vs ZVertex vs #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 399   AttachSPDAcceptance(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NBkgTrackletsVsZVertexVsEta");
 400
 401   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta","Number of charged particles vs ZVertex vs #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 402   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsRecoZVertexVsEta","Number of gen. charged particles vs reco. ZVertex vs gen. #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 403   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",Form("%s vs ZVertex vs #eta;ZVertex;#eta",TrackletsCorrectedName.Data()),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 404   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta");
 405   AttachSPDAcceptance(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsRecoZVertexVsEta");
 406
 407   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsZVertexVsPhi","Number of charged particles vs ZVertex vs #phi;ZVertex;#phi",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),nphibins,phimin,phimax);
 408
 409   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta","Effective number of events vs ZVertex vs #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax()); // Fill 1 unit in each "touched" eta bin per event (represents the eta bins in which each event contributes)
 410   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta");
 411
 412   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Nch",Form("Number of charged particles in |#eta| < %1.1f distribution;N_{ch};N_{events}",fEtaMax),nbinsMult,multMin,multMax);
 413   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"Nch",Form("%s in |#eta| < %1.1f distribution;%s;N_{events}",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax);
 414
 415   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdEta",Form("<dNchd#eta> in |#eta| < %1.1f distribution;dN_{ch}/d#eta;N_{events}",fEtaMax),nbinsMult,multMin,multMax);
 416    CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"dNchdEta",Form("<d%sd#eta> in |#eta| < %1.1f distribution;d%s/d#eta;N_{events}",TrackletsCorrectedAxisName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax);
 417
 418   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsNch",Form("Number of tracklets vs number of generated charged particles in |#eta| < %1.1f;N_{ch};N_{tracklets}",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax); //Response matrix
 419   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"TrackletsVsNch",Form("Number of tracklets vs %s in |#eta| < %1.1f;%s;N_{tracklets}",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 420
 421    CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdetaVsMCdNchdeta",Form("Corrected dN_{ch}/d#eta vs MC dN_{ch}/d#eta in |#eta| < %1.1f;(dN_{ch}/d#eta)_{MC};dN_{ch}/d#eta",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 422
 423   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"V0AMultVsNch",Form("V0A multiplicity vs number of generated charged particles in |#eta| < %1.1f;N_{ch};V0A Mult",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 424   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"V0AMultVsNch",Form("V0A multiplicity vs %s in |#eta| < %1.1f;%s;V0A Mult",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 425
 426   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"V0CMultVsNch",Form("V0C multiplicity vs number of generated charged particles in |#eta| < %1.1f;N_{ch};V0C Mult",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 427   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"V0CMultVsNch",Form("V0C multiplicity vs %s in |#eta| < %1.1f;%s;V0A Mult",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 428
 429 //  CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsTracklets","Number of generated charged particles vs reco tracklets;N_{tracklets};N_{ch}",nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 430
 431   if ( fSPDMeanTrackletsCorrToCompare )
 432   {
 433     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"CorrTrackletsEtaSecVsCorrTrackletsEtaPrim",Form("N_{tr}^{corr} in |#eta|< %1.1f vs N_{tr}^{corr} in |#eta|< %1.1f;N_{tracklets}^{|#eta|<%1.1f};N_{tracklets}^{|#eta|<%1.1f}",fEtaMaxToCompare,fEtaMax,fEtaMax,fEtaMaxToCompare),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 434     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanNchEtaSecVsZVertex",Form("Mean number of corrected tracklets in |#eta| < %1.1f vs Z vertex;Z vertex;<N_{tr}^{corr}>",fEtaMaxToCompare),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 435     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsEtaSecVsZVertex",Form("Mean number of tracklets in |#eta| < %1.1f vs Z vertex;Z vertex;<N_{tr}^{corr}>",fEtaMaxToCompare),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 436     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"TrackletsSecVsZVertexVsEta","Number of tracklets vs ZVertex vs #eta;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 437   }
 438
 439   if ( fSPDMeanTracklets && fSPDOneOverAccxEff )
 440   {
 441     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"dNchdetaComparison2Corrections","dN_{ch}/d#eta comparison with SPD AccxEff and data-driven corrections;dN_{ch}/d#eta;A*N^{corr}_{Tracklets}",151,multMin,150.5,151,multMin,150.5);
 442     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"DispersiondNchdetaComparison2Corrections","dN_{ch}/d#eta dispersion with SPD AccxEff and data-driven corrections;(dN_{ch}/d#eta - A*N^{corr}_{Tracklets}) /dN_{ch}/d#eta;N_{events}",201,-1.005,1.005);
 443     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"CheckMeanNtrCorrVsZVertex","Check for NtrCorr",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 444     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"CheckNtrCorr","Check for NtrCorr distribution;N^{corr}_{tr};N_{events}",nbinsMult,multMin,multMax);
 445   }
 446
 447     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"RelDispersiondNchdetaFromNtrCorrVsdNchdEtaMC",Form("Relative dispersion of dN_{ch}/d#eta from N_{tr}^{corr} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;((dN_{ch}/d#eta)_{gen} - A*N^{corr}_{Tracklets}) /(dN_{ch}/d#eta)_{gen}",fEtaMax),201,-1.005,1.005);
 448     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"DispersiondNchdetaFromNtrCorrVsdNchdEtaMC",Form("Dispersion of dN_{ch}/d#eta from N_{tr}^{corr} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;(dN_{ch}/d#eta)_{gen} - A*N^{corr}_{Tracklets}",fEtaMax),402,-100.5,100.5);
 449     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdetaFromNtrCorrVsdNchdEtaMC",Form("dN_{ch}/d#eta from N_{tr}^{corr} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;(dN_{ch}/d#eta)_{gen};dN_{ch}/d#eta",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 450
 451     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"RelDispersiondNchdetaFromAccEffVsdNchdEtaMC",Form("Relative dispersion of dN_{ch}/d#eta from N_{tr}^{AccxEff} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;((dN_{ch}/d#eta)_{gen} - dN_{ch}/d#eta) /(dN_{ch}/d#eta)_{gen};N_{events}",fEtaMax),201,-1.005,1.005);
 452      CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"DispersiondNchdetaFromAccEffVsdNchdEtaMC",Form("Dispersion of dN_{ch}/d#eta from N_{tr}^{AccxEff} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;(dN_{ch}/d#eta)_{gen} - dN_{ch}/d#eta;N_{events}",fEtaMax),402,-100.5,100.5);
 453     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdetaFromAccEffVsdNchdEtaMC",Form("dN_{ch}/d#eta from N_{tr}^{AccxEff} vs MC dN_{ch}/d#eta in |#eta| < %1.1f;(dN_{ch}/d#eta)_{gen};dN_{ch}/d#eta",fEtaMax),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 454
 455
 456 //   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"CorrTrackletsVsNch","Reco tracklets corrected vs number of generated charged particles;N_{ch};N_{tracklets}^{corr}",nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 457 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"CorrTrackletsVsNch",Form("%s vs number of generated charged particles in |#eta| < %1.1f;N_{ch};%s",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 458
 459   // profile histograms
 460   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta","Mean number of tracklets vs #eta;#eta;<N_{Tracklets}>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0);
 461   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeandNchdEtaVsEta","Mean dN_{ch}/d#eta vs #eta;#eta;<dN_{ch}/d#eta>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0);
 462
 463   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsVsZVertex",Form("Mean number of tracklets in |#eta| < %1.1f vs Z vertex;Z vertex;<N_{Tracklets}>",fEtaMax),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 464
 465   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MeanNchVsEta","Mean number of generated charged particles vs #eta;#eta;<N_{ch}>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0); // Each bin has to be divided by the binwidth to became dNch/dEta (Done in the terminate and stored in MeandNchdEta histo)
 466   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanNchVsEta",Form("Mean %s vs #eta;#eta;<%s>",TrackletsCorrectedName.Data(),TrackletsCorrectedAxisName.Data()),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0); // Each bin has to be divided by the binwidth to became dNch/dEta (Done in the terminate and stored in MeandNchdEta histo)
 467
 468   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MeanNchVsZVertex",Form("Mean number of generated charged particles in |#eta| < %1.1f vs Z vertex;Z vertex;<N_{ch}>",fEtaMax),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 469   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanNchVsZVertex",Form("Mean %s in |#eta| < %1.1f vs Z vertex;Z vertex;<%s>",TrackletsCorrectedName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 470   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeandNchdEtaVsZVertex",Form("Mean d%s/d#eta in |#eta| < %1.1f vs Z vertex;Z vertex;<%s>",TrackletsCorrectedAxisName.Data(),fEtaMax,TrackletsCorrectedAxisName.Data()),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 471
 472
 473   TObjArray* bins = Binning()->CreateBinObjArray("psi","dnchdeta","");
 474   TIter nextBin(bins);
 475   AliAnalysisMuMuBinning::Range* r;
 476
 477   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 478   {
 479     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,Form("EventsIn%s",r->AsString().Data()),Form("Number of %s events in %s bin",triggerClassName,r->AsString().Data()),1,0.,2);
 480   }
 481
 482   delete bins;
 483
 484   TObjArray* binsntr = Binning()->CreateBinObjArray("psi","ntrcorr","");
 485   TIter nextBinNtr(binsntr);
 486   AliAnalysisMuMuBinning::Range* rNtr;
 487
 488   while ( ( rNtr = static_cast<AliAnalysisMuMuBinning::Range*>(nextBinNtr()) ) )
 489   {
 490     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,Form("MeanNchVsZVertex%s",rNtr->AsString().Data()),Form("Mean %s in |#eta| < %1.1f vs Z vertex in bin %s;Z vertex;<%s>",TrackletsCorrectedName.Data(),fEtaMax,rNtr->AsString().Data(),TrackletsCorrectedAxisName.Data()),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 491
 492     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,Form("MeanTrackletsVsZVertex%s",rNtr->AsString().Data()),Form("Mean number of tracklets in |#eta| < %1.1f vs Z vertex in bin %s;Z vertex;<N_{Tracklets}>",fEtaMax,rNtr->AsString().Data()),fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 493   }
 494
 495   delete binsntr;
 496
 497   TObjArray* binsRel = Binning()->CreateBinObjArray("psi","relntrcorr","");
 498   TIter nextBinRel(binsRel);
 499   AliAnalysisMuMuBinning::Range* rRel;
 500
 501   while ( ( rRel = static_cast<AliAnalysisMuMuBinning::Range*>(nextBinRel()) ) )
 502   {
 503     CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,Form("EventsIn%s",rRel->AsString().Data()),Form("Number of %s events in %s bin",triggerClassName,rRel->AsString().Data()),1,0.,2);
 504   }
 505
 506   delete binsRel;
 507
 508 }
 509
 510 //_____________________________________________________________________________
 511 void AliAnalysisMuMuNch::DefineSPDAcceptance()
 512 {
 513   // Defines the functions ( eta = f(z) ) of the edges (right/left) of the inner and outer SPD layers
 514   // R_in = 3.9 cm ; R_out = 7.6 cm
 515
 516   fSPD1LR = new TF1("fSPD1LR",SPDgeomR,-40,40,1);
 517   fSPD1LR->SetParameter(0,3.9);
 518   fSPD1LL = new TF1("fSPD1LL",SPDgeomL,-40,40,1);
 519   fSPD1LL->SetParameter(0,3.9);
 520   fSPD2LR = new TF1("fSPD2LR",SPDgeomR,-40,40,1);
 521   fSPD2LR->SetParameter(0,7.6);
 522   fSPD2LL = new TF1("fSPD2LL",SPDgeomL,-40,40,1);
 523   fSPD2LL->SetParameter(0,7.6);
 524
 525 }
 526
 527 //_____________________________________________________________________________
 528 void AliAnalysisMuMuNch::DefineSPDCorrectionMap(TObjArray* spdCorrectionList)
 529 {
 530   // Defines a TMap of the SPD periods and the corresponding correction
 531   // Each SPD correction in the list must contain in the name the 1st run for which is valid and respect the naming convention:
 532   // (SPDCorrection_1stValidRunNumber_lastValidRunNumber).
 533   // The list must be ordered from smaller to bigger run number
 534   // This is usable with both types of corrections
 535
 536   if ( !fSPDCorrectionMap ) // Creates a TMap for the run->Correction correspondance
 537   {
 538     fSPDCorrectionMap = new TMap;
 539     fSPDCorrectionMap->SetOwnerKeyValue(kTRUE,kTRUE);
 540   }
 541
 542   TIter next(spdCorrectionList);
 543   TH1* SPDCorrection(0x0);
 544   Int_t runRef(0);
 545   Int_t i(0); // SPDCorrection index in the list
 546   while ( (SPDCorrection = static_cast<TH1*>(next())) ) // Checks if the correction list format is ok and if it is ordered
 547   {
 548     if ( static_cast<TH1*>(SPDCorrection)->IsA() != TH2::Class() && static_cast<TH1*>(SPDCorrection)->IsA() != TProfile::Class() )
 549     {
 550       AliFatal("Unrecognized SPD correction Class");
 551     }
 552
 553     TString name = SPDCorrection->GetName();
 554     if ( !name.BeginsWith("SPDCorrection_") ) AliFatal(Form("Incorrect SPD correction at %d format: Objects in list must named as 'SPDCorrection_1stValidRunNumber_lastValidRunNumber'",i));
 555
 556     name.Remove(0,name.First("_") + 1);
 557     TString nameLast = name;
 558     name.Remove(name.First("_"),name.Length());
 559     nameLast.Remove(0,nameLast.First("_") + 1);
 560     if ( !name.IsDigit() ) AliFatal(Form("Incorrect SPD correction at %d format: Impossible to retrieve first valid run number",i));
 561     if ( !nameLast.IsDigit() ) AliFatal(Form("Incorrect SPD correction at %d format: Impossible to retrieve last valid run number",i));
 562
 563     Int_t runLow = name.Atoi();
 564     Int_t runHigh = nameLast.Atoi();
 565     if ( runHigh < runLow ) AliFatal(Form("SPD correction at %d validity range not valid",i));
 566     if ( runLow <= runRef ) AliFatal("SPD correction list not in correct order. Should be ordered from smaller to bigger run number");
 567     else runRef = runHigh;
 568
 569     for ( Int_t j = runLow ; j <= runHigh ; j++ )
 570     {
 571       fSPDCorrectionMap->Add(new TObjString(Form("%d",j)),new TObjString(Form("%d",i))); // The mapping is done between the runs and the index of the correction in the spdCorrection List
 572     }
 573
 574     i++;
 575   }
 576
 577   fSPDCorrectionList = static_cast<TObjArray*>(spdCorrectionList->Clone());
 578 }
 579
 580 //_____________________________________________________________________________
 581 void AliAnalysisMuMuNch::AddHisto(const char* eventSelection,
 582                                   const char* triggerClassName,
 583                                   const char* centrality,
 584                                   const char* histoname,
 585                                   Double_t z,
 586                                   TH1* h,Bool_t isMC)
 587 {
 588   // Adds the content of a 1D histo to a 2D histo a the z position
 589
 590   Int_t zbin = fZAxis->FindBin(z);
 591   TH2F* h2;
 592   if (isMC) h2 = static_cast<TH2F*>(MCHisto(eventSelection,triggerClassName,centrality,histoname));
 593   else h2 = static_cast<TH2F*>(Histo(eventSelection,triggerClassName,centrality,histoname));
 594
 595   for ( Int_t i = 1; i <= h->GetXaxis()->GetNbins(); ++i )
 596   {
 597     Double_t content = h2->GetCellContent(zbin,i);
 598
 599     if ( h->GetBinContent(i) >  0 )
 600     {
 601       h2->SetCellContent(zbin,i,content + h->GetBinContent(i));
 602     }
 603   }
 604
 605   h2->SetEntries(h2->GetSumOfWeights());
 606 }
 607
 608
 609 //_____________________________________________________________________________
 610 void AliAnalysisMuMuNch::AttachSPDAcceptance(UInt_t dataType,
 611                                              const char* eventSelection,
 612                                              const char* triggerClassName,
 613                                              const char* centrality,const char* histoname)
 614 {
 615   // Attachs AccxEff curves to the histogram
 616
 617   if ( dataType & kHistoForData )
 618   {
 619     if( !Histo(eventSelection,triggerClassName,centrality,histoname) )
 620     {
 621       AliError(Form("ERROR: SPD Acceptance curves attach failed. Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
 622       return;
 623     }
 624
 625     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LR->Clone());
 626     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LL->Clone());
 627     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LR->Clone());
 628     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LL->Clone());
 629   }
 630   if ( (dataType & kHistoForMCInput) && HasMC() )
 631   {
 632     if( !MCHisto(eventSelection,triggerClassName,centrality,histoname) )
 633     {
 634       AliError(Form("ERROR: SPD Acceptance curves attach failed. MC Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
 635       return;
 636     }
 637
 638     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LR->Clone());
 639     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LL->Clone());
 640     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LR->Clone());
 641     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LL->Clone());
 642
 643   }
 644
 645 }
 646
 647
 648 //_____________________________________________________________________________
 649 void AliAnalysisMuMuNch::FillHistosForEvent(const char* eventSelection,
 650                                             const char* triggerClassName,
 651                                             const char* centrality)
 652 {
 653   // Fills the data (or reco if simu) multiplicity histos
 654
 655   if ( IsHistogrammingDisabled() ) return;
 656
 657   if ( fResolution ) return; //When computing resolutions we skip this method
 658
 659 //  if ( !AliAnalysisMuonUtility::IsAODEvent(Event()) )
 660 //  {
 661 //    AliError("Don't know how to deal with ESDs...");
 662 //    return;
 663 //  }
 664
 665 //  if ( HasMC() && !fSPDOneOverAccxEff && !fSPDMeanTracklets ) // We have MC but no correction (SPD correction computation mode)so we skip the method
 666 //  {
 667 //    return;
 668 //  }
 669
 670 //  AliAODEvent* aod = static_cast<AliAODEvent*>(Event());
 671 //
 672 //  AliVVertex* vertex = aod->GetPrimaryVertexSPD();
 673
 674   AliVVertex* vertex = AliAnalysisMuonUtility::GetVertexSPD(Event());
 675
 676   TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
 677
 678   if (!nchList || nchList->IsEmpty() ) // Empty NCH means that there is no SPD vertex ( see SetEvent() ) when runing on data.
 679   {
 680     AliError("Empty Nch list in event");
 681     return;
 682   }
 683 //  nchList->Print();
 684   Double_t SPDZv;
 685
 686   if ( !vertex || vertex->GetZ() == 0.0 ) // Running in Simu the spdZ == 0 means no SPD info. In data, avoid breaks in events w/o SPD vertex
 687   {
 688     SPDZv = -40.;
 689   }
 690   else SPDZv = vertex->GetZ();
 691
 692   TH1* hSPDcorrectionVsEta = static_cast<TH1*>(nchList->FindObject("SPDcorrectionVsEta"));
 693   TH1* hNTrackletVsEta = static_cast<TH1*>(nchList->FindObject("NTrackletVsEta"));
 694   TH1* hNTrackletVsPhi = static_cast<TH1*>(nchList->FindObject("NTrackletVsPhi"));
 695
 696   TH1* hNTrackletSecVsEta(0x0);
 697   if ( fSPDMeanTrackletsCorrToCompare ) // In case we have a secondary eta range we clone the tracklets vs eta histo
 698   {
 699     hNTrackletSecVsEta = static_cast<TH1*>(hNTrackletVsEta->Clone());
 700     hNTrackletSecVsEta->Reset();
 701   }
 702
 703   TH1* hNchVsEta = static_cast<TH1*>(hNTrackletVsEta->Clone("NchVsEta"));
 704
 705   TProfile* hMeanTrackletsVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta"));
 706   TProfile* hMeanNchVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
 707   TProfile* hMeandNchdEtaVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeandNchdEtaVsEta"));
 708
 709   TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(Histo(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
 710
 711   Int_t nBins(0);
 712
 713   Int_t nTracklets[2] = {0,0}; // {nTracklets_Eta_fetaRange,nTracklets_Eta_fetaRangeComp} The first component is the number of tracklets in the primary eta range and the second is  the number of tracklets in the secondary eta range
 714   Double_t nch[2] = {0.,0.}; // {nCorrTracklets_Eta_fetaRange,nCorrTracklets_Eta_fetaRangeComp} The first component is the corrected(with the primary correction) number of tracklets in the primary eta range and the second is  the number of corrected(with the secondary correction) tracklets in the secondary eta range
 715
 716   Int_t etaBinMin[2] = {fEtaAxis->FindBin(fetaRange[0]),fEtaAxis->FindBin(fEtaMinToCompare)}; // Primary and secondary (if any) eta ranges
 717   Int_t etaBinMax[2] = {fEtaAxis->FindBin(fetaRange[1]),fEtaAxis->FindBin(fEtaMaxToCompare)-1}; // Note that the binMax is the bin-1 because the upper extreme of a bin belongs to the next bin (i.e. the bin containing -0.5 has as center -0.45, the bin is [-0.5,-0.4) but the bin containing 0.5 has as center 1.05, the bin is [0.5,1.1) so we have to take the previous bin [0.4,0.5)) This is already taken into account in the GetEtaRangeSPD method
 718
 719   for (Int_t j = etaBinMin[0] ; j <= etaBinMax[0] ; j++) // Loop over eta bins.
 720   {
 721     Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
 722
 723     Double_t eta = fEtaAxis->GetBinCenter(j);
 724
 725     if ( correction < 0 ) continue;
 726     else if ( correction == 0.0 ) // No tracklets found in this eta bin.
 727     {
 728       correction = GetSPDCorrection(SPDZv,eta); //FIXME: Here eta is the bincenter not the exact one, will be a problem if the SPD AccxEff binning is not the same as the fEtaAxis one
 729
 730       if ( correction == 0. || correction > 2.5) continue; // We need to know if the eta bin is in a region we want to take into account to count or not the zero
 731     }
 732
 733     Int_t ntr = hNTrackletVsEta->GetBinContent(j); // Tracklets in eta bin
 734     nTracklets[0] += ntr;
 735     if ( j >= etaBinMin[1] && j <= etaBinMax[1] )
 736     {
 737       nTracklets[1] += ntr; // number of tracklets in the secondary eta range
 738       hNTrackletSecVsEta->SetBinContent(j,ntr);
 739     }
 740
 741     hMeanTrackletsVsEta->Fill(eta,ntr); // Fill the number of tracklets of each eta bin in the profile
 742
 743     if ( fSPDOneOverAccxEff )
 744     {
 745       nch[0] += ntr * correction; // Number of charged particles (SPD AccxEff corrected tracklets)
 746       hMeanNchVsEta->Fill(eta,ntr*correction);
 747       hMeandNchdEtaVsEta->Fill(eta,ntr*correction/fEtaAxis->GetBinWidth(5));
 748       hNchVsEta->SetBinContent(j,hNchVsEta->GetBinContent(j) * correction ); // Fill the number of charged particles of each eta bin in the profile
 749     }
 750
 751     ++nBins; // We sum up the number of bins entering in the computation
 752
 753     hEventsVsZVertexVsEta->Fill(SPDZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
 754   }
 755
 756   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi",SPDZv,hNTrackletVsPhi);
 757   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta",SPDZv,hNTrackletVsEta);
 758
 759   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",SPDZv,hNchVsEta);
 760
 761   Histo(eventSelection,triggerClassName,centrality,"Tracklets")->Fill(nTracklets[0]);
 762   Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsVsZVertex")->Fill(SPDZv,nTracklets[0]);
 763
 764
 765   if ( fSPDMeanTracklets && !fSPDOneOverAccxEff ) // Just correct by SPDMeanTracklets method if there is no SPD AccxEff correction
 766   {
 767     Double_t SPDr = GetTrackletsMeanCorrection(SPDZv,nTracklets[0]); // Get 'mean correction' for the zvtx
 768
 769     if ( SPDr < -999.) nch[0] = -1;
 770     else nch[0] = nTracklets[0] + SPDr; // In case of 'mean correction' nch has not be filled in the eta bins loop
 771
 772     if ( fSPDMeanTrackletsCorrToCompare ) // Comparison of corrected tracklets in the primary and secondary eta ranges
 773     {
 774       AddHisto(eventSelection,triggerClassName,centrality,"TrackletsSecVsZVertexVsEta",SPDZv,hNTrackletSecVsEta);
 775
 776       Double_t SPDrEtaComp = GetTrackletsMeanCorrection(SPDZv,nTracklets[1],kTRUE); // Get secondary 'mean correction' for the zvtx
 777
 778       if ( SPDrEtaComp < -999.) nch[1] = -1;
 779       else nch[1] = nTracklets[1] + SPDrEtaComp; // In case of 'mean correction' nch has not be filled in the eta bins loop
 780
 781       Histo(eventSelection,triggerClassName,centrality,"CorrTrackletsEtaSecVsCorrTrackletsEtaPrim")->Fill(nch[0],nch[1]);
 782       Histo(eventSelection,triggerClassName,centrality,"MeanNchEtaSecVsZVertex")->Fill(SPDZv,nch[1]); // Control plot to check if the secondary correction is applied correctly
 783       Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsEtaSecVsZVertex")->Fill(SPDZv,nTracklets[1]);
 784     }
 785
 786   }
 787   else if ( fSPDMeanTracklets && fSPDOneOverAccxEff ) // Here we compare the Corrected tracklets with the two correction methods
 788   {
 789     Int_t i(-1);
 790     Bool_t parFound1(kFALSE),parFound2(kFALSE);
 791     Double_t NtrCorr(0.),dNchdeta(0.);
 792     while ( i < nchList->GetEntries() - 1 && !(parFound1 && parFound2) )
 793     {
 794       i++;
 795       while ( nchList->At(i)->IsA() != TParameter<Double_t>::Class() && i < nchList->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
 796       {
 797         i++;
 798       }
 799
 800       TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(nchList->At(i));
 801
 802       if ( TString(p->GetName()).Contains("NtrCorr") )
 803       {
 804         parFound1 = kTRUE;
 805         NtrCorr = p->GetVal();
 806       }
 807       else if ( TString(p->GetName()).Contains("MeandNchdEta") )
 808       {
 809         parFound2 = kTRUE;
 810         dNchdeta = p->GetVal();
 811       }
 812     }
 813
 814     if ( parFound1 && parFound2 )
 815     {
 816 //      Double_t dNchdetaPubli = 17.35; //FIXME: hardcoded (pPb value)
 817       Double_t ctToNch = 1.11; //FIXME: hardcoded (value for Nch vs NtrCorr(eta<0.5) in pPb)
 818
 819       Histo(eventSelection,triggerClassName,centrality,"dNchdetaComparison2Corrections")->Fill(dNchdeta,ctToNch*NtrCorr/(2*fEtaMax));
 820       Histo(eventSelection,triggerClassName,centrality,"CheckMeanNtrCorrVsZVertex")->Fill(SPDZv,NtrCorr);
 821       if ( dNchdeta !=0 )
 822       {
 823         Histo(eventSelection,triggerClassName,centrality,"DispersiondNchdetaComparison2Corrections")->Fill((dNchdeta - ctToNch*NtrCorr/(2*fEtaMax)) / dNchdeta);
 824       }
 825
 826       Histo(eventSelection,triggerClassName,centrality,"CheckNtrCorr")->Fill(NtrCorr);
 827
 828     }
 829   }
 830
 831
 832   Histo(eventSelection,triggerClassName,centrality,"TrackletsVsNch")->Fill(nch[0],nTracklets[0]);
 833   Histo(eventSelection,triggerClassName,centrality,"Nch")->Fill(nch[0]);
 834   Histo(eventSelection,triggerClassName,centrality,"MeanNchVsZVertex")->Fill(SPDZv,nch[0]);
 835
 836
 837   Double_t V0AMult = 0.;
 838   Double_t V0CMult = 0.;
 839
 840 //  AliVVZERO* vzero = aod->GetVZEROData();
 841   AliVVZERO* vzero = Event()->GetVZEROData();
 842   if (vzero)
 843   {
 844     Double_t multV0A = vzero->GetMTotV0A();
 845     V0AMult = AliESDUtils::GetCorrV0A(multV0A,SPDZv);
 846     Double_t multV0C = vzero->GetMTotV0C();
 847     V0CMult = AliESDUtils::GetCorrV0C(multV0C,SPDZv);
 848
 849     Histo(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nch[0]);
 850     Histo(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nch[0]);
 851
 852   }
 853
 854
 855   delete hNchVsEta;
 856
 857   // Mean dNch/dEta computation. In case the correction is the 'mean' one, this is not really dNch/deta, we have to multiply it by a factor <Nch>/<Ntrkls_corr> (we get <Nch> from a MC)
 858   Double_t meandNchdEta(0.);
 859
 860   if ( nBins >  0 )
 861   {
 862     if ( fSPDOneOverAccxEff ) meandNchdEta = nch[0] / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
 863     else if ( fSPDMeanTracklets ) meandNchdEta = nch[0] / (2.*fEtaMax); //fEtaAxis->GetBinWidth(5);
 864   }
 865
 866   Histo(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
 867   Histo(eventSelection,triggerClassName,centrality,"MeandNchdEtaVsZVertex")->Fill(SPDZv,meandNchdEta);
 868
 869
 870   //_____________These were tests //FIXME: Check if this tests are still neccesary_____________
 871   TObjArray* bins = Binning()->CreateBinObjArray("psi","dnchdeta","");
 872   TIter nextBin(bins);
 873   AliAnalysisMuMuBinning::Range* r;
 874
 875   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 876   {
 877     Bool_t ok(kFALSE);
 878     if ( r->Quantity() == "DNCHDETA" )
 879     {
 880       ok = r->IsInRange(meandNchdEta);
 881     }
 882
 883     if ( ok )
 884     {
 885       Histo(eventSelection,triggerClassName,centrality,Form("EventsIn%s",r->AsString().Data()))->Fill(1.);
 886     }
 887   }
 888
 889   delete bins;
 890
 891
 892   TObjArray* binsNtrRaw = Binning()->CreateBinObjArray("psi","ntrraw","");
 893   TIter nextBinNtrRaw(binsNtrRaw);
 894   AliAnalysisMuMuBinning::Range* rNtrRaw;
 895
 896   while ( ( rNtrRaw = static_cast<AliAnalysisMuMuBinning::Range*>(nextBinNtrRaw()) ) )
 897   {
 898     Bool_t ok(kFALSE);
 899     if ( rNtrRaw->Quantity() == "NTRRAW" )
 900     {
 901       ok = rNtrRaw->IsInRange(nTracklets[0]);
 902     }
 903
 904     if ( ok )
 905     {
 906       Histo(eventSelection,triggerClassName,centrality,Form("MeanTrackletsVsZVertex%s",rNtrRaw->AsString().Data()))->Fill(SPDZv,nTracklets[0]);
 907     }
 908   }
 909
 910   delete binsNtrRaw;
 911
 912
 913   TObjArray* binsNtr = Binning()->CreateBinObjArray("psi","ntrcorr","");
 914   TIter nextBinNtr(binsNtr);
 915   AliAnalysisMuMuBinning::Range* rNtr;
 916
 917   while ( ( rNtr = static_cast<AliAnalysisMuMuBinning::Range*>(nextBinNtr()) ) )
 918   {
 919     Bool_t ok(kFALSE);
 920     if ( rNtr->Quantity() == "NTRCORR" )
 921     {
 922       ok = rNtr->IsInRange(nch[0]);
 923     }
 924
 925     if ( ok )
 926     {
 927       Histo(eventSelection,triggerClassName,centrality,Form("MeanNchVsZVertex%s",rNtr->AsString().Data()))->Fill(SPDZv,nch[0]);
 928     }
 929   }
 930
 931   delete binsNtr;
 932
 933
 934   TObjArray* binsRel = Binning()->CreateBinObjArray("psi","relntrcorr","");
 935   TIter nextBinRel(binsRel);
 936   AliAnalysisMuMuBinning::Range* rRel;
 937   while ( ( rRel = static_cast<AliAnalysisMuMuBinning::Range*>(nextBinRel()) ) )
 938   {
 939     Bool_t ok(kFALSE);
 940     if ( rRel->Quantity() == "RELNTRCORR" )
 941     {
 942       ok = rRel->IsInRange(nch[0]/fMeanTrRef);
 943     }
 944
 945     if ( ok )
 946     {
 947       Histo(eventSelection,triggerClassName,centrality,Form("EventsIn%s",rRel->AsString().Data()))->Fill(1.);
 948     }
 949   }
 950
 951   delete binsRel;
 952   //_________________________________________________________________________________________
 953
 954 }
 955
 956 //_____________________________________________________________________________
 957 void AliAnalysisMuMuNch::FillHistosForMCEvent(const char* eventSelection,const char* triggerClassName,const char* centrality)
 958 {
 959   /// Fill input MC multiplicity histos
 960
 961   if ( IsHistogrammingDisabled() ) return;
 962
 963   TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
 964
 965   if (!nchList || nchList->IsEmpty())
 966   {
 967     AliError("Empty Nch list in event");
 968     return;
 969   }
 970
 971   Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent()); // Definition of MC generated z vertex
 972 //  AliVVertex* vertex = static_cast<AliAODEvent*>(Event())->GetPrimaryVertexSPD();
 973   AliVVertex* vertex = AliAnalysisMuonUtility::GetVertexSPD(Event());
 974
 975   Double_t SPDZv(0.);
 976   TParameter<Double_t>* p(0x0);
 977
 978   //____Resolution Histos___
 979   if ( !fSPDOneOverAccxEff && !fSPDMeanTracklets && fResolution )
 980   {
 981     Int_t nContributors(0);
 982     if ( vertex )
 983     {
 984       nContributors  = vertex->GetNContributors();
 985       SPDZv = vertex->GetZ();
 986     }
 987     MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsnC")->Fill(nContributors,SPDZv - MCZv);
 988     MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsMCz")->Fill(MCZv,SPDZv - MCZv);
 989
 990     Double_t EtaReco(0.),EtaMC(0.),PhiReco(0.),PhiMC(0.);
 991     Int_t i(-1),labelEtaReco(-1),labelEtaMC(-1),labelPhiReco(-1),labelPhiMC(-1);
 992
 993     while ( i < nchList->GetEntries() - 1 )
 994     {
 995       i++;
 996       while ( nchList->At(i)->IsA() != TParameter<Double_t>::Class() ) // In case there is a diferent object, just to skip it
 997       {
 998         i++;
 999       }
1000
1001       p = static_cast<TParameter<Double_t>*>(nchList->At(i));
1002
1003
1004          //Eta Resolution
1005       if ( TString(p->GetName()).Contains("EtaReco") ) // We take the reco eta
1006       {
1007         sscanf(p->GetName(),"EtaReco%d",&labelEtaReco);
1008         EtaReco = p->GetVal();
1009         MCHisto(eventSelection,triggerClassName,centrality,"Eta")->Fill(EtaReco);
1010       }
1011       else if ( TString(p->GetName()).Contains("EtaMC") ) // We take the generated eta
1012       {
1013         sscanf(p->GetName(),"EtaMC%d",&labelEtaMC);
1014         EtaMC = p->GetVal();
1015         MCHisto(eventSelection,triggerClassName,centrality,"MCEta")->Fill(EtaMC);
1016       }
1017       if ( labelEtaReco > 0 && labelEtaReco == labelEtaMC ) // To be sure we compute the difference for the same particle
1018       {
1019         labelEtaReco = -1; // Restart of the label value to avoid double count the eta difference when computing the phi one
1020         Double_t EtaDif = EtaReco - EtaMC;
1021         MCHisto(eventSelection,triggerClassName,centrality,"EtaRes")->Fill(EtaDif);
1022         MCHisto(eventSelection,triggerClassName,centrality,"EtaResVsZ")->Fill(MCZv,EtaDif);
1023         MCHisto(eventSelection,triggerClassName,centrality,"EtaResVsnC")->Fill(nContributors,EtaDif);
1024       }
1025
1026          //Phi Resolution
1027       else if ( TString(p->GetName()).Contains("PhiReco") ) // We take the reco phi
1028       {
1029         sscanf(p->GetName(),"PhiReco%d",&labelPhiReco);
1030         PhiReco = p->GetVal();
1031         MCHisto(eventSelection,triggerClassName,centrality,"Phi")->Fill(PhiReco);
1032       }
1033       else if ( TString(p->GetName()).Contains("PhiMC") ) // We take the generated phi
1034       {
1035         sscanf(p->GetName(),"PhiMC%d",&labelPhiMC);
1036         PhiMC = p->GetVal();
1037         MCHisto(eventSelection,triggerClassName,centrality,"MCPhi")->Fill(PhiMC);
1038       }
1039
1040       if ( labelPhiReco > 0 && labelPhiReco == labelPhiMC ) // To be sure we compute the difference for the same particle
1041       {
1042         labelPhiReco = -1; // Restart of the label value to avoid double count the phi difference when computing the eta one
1043         Double_t PhiDif = PhiReco - PhiMC;
1044         MCHisto(eventSelection,triggerClassName,centrality,"PhiRes")->Fill(PhiDif);
1045
1046         //___With the following algorithm we refer the differences to the interval [-Pi/2,Pi/2]
1047         if ( PhiDif < -TMath::PiOver2() && PhiDif > -TMath::Pi() )
1048         {
1049           PhiDif = -TMath::Pi() - PhiDif;
1050         }
1051         else if ( PhiDif < -TMath::Pi() )
1052         {
1053           PhiDif = -2.*TMath::Pi() - PhiDif;
1054           if ( PhiDif < -TMath::PiOver2() )
1055           {
1056             PhiDif = -TMath::Pi() - PhiDif;
1057           }
1058         }
1059
1060         else if ( PhiDif > TMath::PiOver2() && PhiDif < TMath::Pi() )
1061         {
1062           PhiDif = TMath::Pi() - PhiDif;
1063         }
1064         else if ( PhiDif > TMath::Pi() )
1065         {
1066           PhiDif = 2.*TMath::Pi() - PhiDif;
1067           if ( PhiDif > TMath::PiOver2() )
1068           {
1069             PhiDif = TMath::Pi() - PhiDif;
1070           }
1071         }
1072         //___
1073
1074         MCHisto(eventSelection,triggerClassName,centrality,"PhiResVsZ")->Fill(MCZv,PhiDif);
1075         MCHisto(eventSelection,triggerClassName,centrality,"PhiResShifted")->Fill(PhiDif);
1076         MCHisto(eventSelection,triggerClassName,centrality,"PhiResVsnC")->Fill(nContributors,PhiDif);
1077       }
1078     }
1079
1080     return; // When computing resolutions we skip the rest of the method
1081   }
1082   //_____
1083
1084   //___Multiplicity histos (just when not computing resolutions)___
1085   TH1* hSPDcorrectionVsEta = static_cast<TH1*>(nchList->FindObject("MCSPDcorrectionVsEta"));
1086   TH1* hNBkgTrackletsVSEta = static_cast<TH1*>(nchList->FindObject("NBkgTrackletsVSEta"));
1087   TH1* hNchVsPhi = static_cast<TH1*>(nchList->FindObject("MCNchVsPhi"));
1088   TH1* hNchVsEta = static_cast<TH1*>(nchList->FindObject("MCNchVsEta"));
1089   TH1* hNTrackletVsEta = static_cast<TH1*>(nchList->FindObject("MCNTrackletVsEta"));
1090   TH1* hNTrackletVsPhi = static_cast<TH1*>(nchList->FindObject("MCNTrackletVsPhi"));
1091
1092   TProfile* hMeanNchVsEta = static_cast<TProfile*>(MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
1093
1094   TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(MCHisto(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
1095
1096   Int_t nBins(0);
1097
1098   Double_t nchSum(0.),nTracklets(0.);
1099   for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) //Loop over all eta bins
1100   {
1101     Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
1102
1103     if ( correction < 0 ) continue; // We count just the particles in the SPD acceptance. (MCSPDcorrectionVsEta is filled to -1 in SetEvent() for those bins out of range)
1104
1105     nTracklets += hNTrackletVsEta->GetBinContent(j); // Reco tracklets
1106
1107     ++nBins; // We sum up the number of bins entering in the computation
1108
1109     Double_t eta = fEtaAxis->GetBinCenter(j);
1110     Double_t nch = hNchVsEta->GetBinContent(j); // Generated particles
1111
1112     nchSum += nch;
1113
1114     hMeanNchVsEta->Fill(eta,nch); // Fill the number of charged particles of each eta bin in the profile
1115
1116     hEventsVsZVertexVsEta->Fill(MCZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
1117   }
1118
1119   MCHisto(eventSelection,triggerClassName,centrality,"Tracklets")->Fill(nTracklets); // Note that these are NOT the same tracklets as in the FillHistosForEvent() since here the SPD "dead" zones (the ones where correction > threshold) are not rejected
1120
1121   if ( vertex )
1122   {
1123     SPDZv = vertex->GetZ();
1124   }
1125
1126   Bool_t isMChisto = kTRUE; // Used to get the MC histos in the Add method
1127
1128   AddHisto(eventSelection,triggerClassName,centrality,"NBkgTrackletsVsZVertexVsEta",SPDZv,hNBkgTrackletsVSEta,isMChisto);
1129   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsPhi",MCZv,hNchVsPhi,isMChisto);
1130   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",MCZv,hNchVsEta,isMChisto);
1131   AddHisto(eventSelection,triggerClassName,centrality,"NchVsRecoZVertexVsEta",SPDZv,hNchVsEta,isMChisto);
1132   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta",SPDZv,hNTrackletVsEta,isMChisto);
1133   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi",SPDZv,hNTrackletVsPhi,isMChisto);
1134
1135   MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsZVertex")->Fill(MCZv,nchSum);
1136   MCHisto(eventSelection,triggerClassName,centrality,"Nch")->Fill(nchSum);
1137
1138
1139   // Mean dNch/dEta computation
1140   Double_t meandNchdEta(0.);
1141   if ( nBins >  0 )
1142   {
1143     meandNchdEta = nchSum / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
1144   }
1145
1146   MCHisto(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
1147
1148
1149   Int_t i(-1);
1150   p = 0x0;
1151   Double_t nTrCorr(0.);
1152   Double_t dNchdetaReco(0.);
1153   Double_t ctToNch = 1.11; //FIXME: hardcoded (value for Nch vs NtrCorr(eta<0.5) in pPb)
1154   while ( i < nchList->GetEntries() - 1 )
1155   {
1156     i++;
1157     while ( nchList->At(i)->IsA() != TParameter<Double_t>::Class() && i < nchList->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
1158     {
1159       i++;
1160     }
1161
1162     p = static_cast<TParameter<Double_t>*>(nchList->At(i));
1163
1164     if ( ( TString(p->GetName()).Contains("NtrCorr") || TString(p->GetName()).BeginsWith("Nch") ) && ( (fSPDMeanTracklets && !fSPDOneOverAccxEff) || (!fSPDMeanTracklets && fSPDOneOverAccxEff) ))
1165     {
1166       MCHisto(eventSelection,triggerClassName,centrality,"CorrTrackletsVsNch")->Fill(nchSum,p->GetVal());
1167     }
1168     else if ( TString(p->GetName()).Contains("NtrCorr") )
1169     {
1170       nTrCorr = p->GetVal();
1171       MCHisto(eventSelection,triggerClassName,centrality,"dNchdetaFromNtrCorrVsdNchdEtaMC")->Fill(meandNchdEta,ctToNch*nTrCorr/(2*fEtaMax));
1172       MCHisto(eventSelection,triggerClassName,centrality,"RelDispersiondNchdetaFromNtrCorrVsdNchdEtaMC")->Fill((meandNchdEta - ctToNch*nTrCorr/(2*fEtaMax)) / meandNchdEta);
1173       MCHisto(eventSelection,triggerClassName,centrality,"DispersiondNchdetaFromNtrCorrVsdNchdEtaMC")->Fill((meandNchdEta - ctToNch*nTrCorr/(2*fEtaMax)));
1174     }
1175     else if ( TString(p->GetName()).Contains("Ntr") && !TString(p->GetName()).Contains("SPDOk") && !TString(p->GetName()).Contains("Corr"))
1176     {
1177       MCHisto(eventSelection,triggerClassName,centrality,"TrackletsVsNch")->Fill(nchSum,p->GetVal());
1178     }
1179     else if ( TString(p->GetName()).Contains("MeandNchdEta") )
1180     {
1181       dNchdetaReco = p->GetVal();
1182       MCHisto(eventSelection,triggerClassName,centrality,"dNchdetaVsMCdNchdeta")->Fill(meandNchdEta,dNchdetaReco);
1183       MCHisto(eventSelection,triggerClassName,centrality,"dNchdetaFromAccEffVsdNchdEtaMC")->Fill(meandNchdEta,dNchdetaReco);
1184       MCHisto(eventSelection,triggerClassName,centrality,"RelDispersiondNchdetaFromAccEffVsdNchdEtaMC")->Fill((meandNchdEta - dNchdetaReco) / meandNchdEta);
1185       MCHisto(eventSelection,triggerClassName,centrality,"DispersiondNchdetaFromAccEffVsdNchdEtaMC")->Fill((meandNchdEta - dNchdetaReco));
1186     }
1187   }
1188
1189   Double_t V0AMult = 0.;
1190   Double_t V0CMult = 0.;
1191
1192   AliVVZERO* vzero = Event()->GetVZEROData();
1193   if (vzero)
1194   {
1195     Double_t multV0A = vzero->GetMTotV0A();
1196     V0AMult = AliESDUtils::GetCorrV0A(multV0A,MCZv);
1197     Double_t multV0C = vzero->GetMTotV0C();
1198     V0CMult = AliESDUtils::GetCorrV0C(multV0C,MCZv);
1199
1200     MCHisto(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nchSum);
1201     MCHisto(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nchSum);
1202
1203   }
1204   //____
1205 }
1206
1207
1208 //_____________________________________________________________________________
1209 Bool_t AliAnalysisMuMuNch::GetEtaRangeSPD(Double_t spdZVertex, Double_t etaRange[])
1210 {
1211   // Fill the SPD eta range for a given z vertex and returns whether the range is valid or not.
1212
1213   if ( fSPDMeanTracklets || (!fSPDMeanTracklets && !fSPDOneOverAccxEff) )
1214   {
1215     // If we are using the mean correction or we want the raw distributions we should not apply the algorithm to restric the eta range, in order to avoid the effect of the eta-Zv bins on the distributions ( the algorithm makes that when an eta bin is partially out of the SPD acceptance we remove it. This makes that the mean raw tracklets distribution is not smooth for Zv in which the eta range is partially out of the SPD )
1216
1217     //FIXME: Implement something to warning or reject events for which the selected eta range is completely out of the SPD acceptance
1218     etaRange[0] = fEtaMin;
1219     etaRange[1] = fEtaMax - 1E-6;
1220
1221     return kTRUE;
1222   }
1223
1224   Double_t etaMax(fEtaMax),etaMin(fEtaMin);
1225
1226   Double_t vf1LR = fSPD1LR->Eval(spdZVertex); //Eta values for the z vertex over the SPD acceptance curves
1227   Double_t vf2LR = fSPD2LR->Eval(spdZVertex);
1228   Double_t vf1LL = fSPD1LL->Eval(spdZVertex);
1229   Double_t vf2LL = fSPD2LL->Eval(spdZVertex);
1230
1231   //____We start by asigning the eta range as the eta values in the SPD acceptance curve
1232   if ( spdZVertex < 0. )
1233   {
1234     etaMax = vf2LR;
1235     etaMin = TMath::Max(vf1LL,vf2LL);
1236   }
1237   else
1238   {
1239     etaMax = TMath::Min(vf1LR,vf2LR);
1240     etaMin = vf2LL;
1241   }
1242   //____
1243
1244
1245   //____Algorithm to avoid taking bins which are crossed by the SPD acceptance curves
1246   Int_t binYMin = fEtaAxis->FindBin(etaMin); // Find the corresponding bins for eta max & min and z vertex
1247   Int_t binYMax = fEtaAxis->FindBin(etaMax);
1248   Int_t binX = fZAxis->FindBin(spdZVertex);
1249
1250   // Define the values for the relevant edges of the eta and z bins
1251   Double_t upEdge = fEtaAxis->GetBinUpEdge(binYMax);  // up edge of the top eta bin
1252   Double_t lowEdge = fEtaAxis->GetBinLowEdge(binYMin); // low edge of the bottom eta bin
1253   Double_t leftEdge = fZAxis->GetBinLowEdge(binX); // left edge of the z bin
1254   Double_t rightEdge = fZAxis->GetBinUpEdge(binX); // right edge of the z bin
1255
1256   Double_t etaMaxTemp(0.),etaMinTemp(0.);
1257   if ( spdZVertex < 0. )
1258   {
1259     etaMaxTemp = fSPD2LR->Eval(rightEdge); // Define the temporary eta max as the value of the curve int the righ edge of the bin
1260     etaMinTemp = TMath::Max(fSPD1LL->Eval(leftEdge),fSPD2LL->Eval(leftEdge));
1261   }
1262   else
1263   {
1264     etaMaxTemp = TMath::Min(fSPD1LR->Eval(rightEdge),fSPD2LR->Eval(rightEdge));
1265     etaMinTemp = fSPD2LL->Eval(leftEdge);
1266   }
1267
1268   while ( upEdge > etaMaxTemp ) //We take eta max as the up edge of the 1st bin which is inside the SPD acceptance but not crossed by the curve
1269   {
1270     binYMax = binYMax - 1; // Since the up edge of the current bin is bigger than the SPD acc curve we move 1 bin below
1271     upEdge = fEtaAxis->GetBinUpEdge(binYMax); // Take the up edge of the new bin
1272     etaMax = upEdge - 1E-6; // We substract 1E-6 cause the up edge of the a bin belongs to the bin+1
1273
1274   }
1275
1276   //We take eta min as the low edge of the 1st bin which is inside the SPD acceptance but not crossed by the curve
1277   while ( lowEdge < etaMinTemp )
1278   {
1279     binYMin = binYMin + 1; // Since the low edge of the current bin is smaller than the SPD acc curve we move 1 bin above
1280     lowEdge = fEtaAxis->GetBinLowEdge(binYMin); // Take the low edge of the new bin
1281     etaMin = lowEdge + 1E-6;
1282   }
1283   //____
1284
1285   // In case the eta range we want (given in the constructor) is smaller than the one found by the algorithm (max range) we redefine the values
1286   if ( etaMin < fEtaMin ) etaMin = fEtaMin + 1E-6;
1287   if ( etaMax > fEtaMax ) etaMax = fEtaMax - 1E-6;
1288
1289   if ( etaMin > etaMax ) // If this happens, the event is not able to cover the desired eta range, so we set an invalid range
1290   {
1291     etaRange[0] = 999.;
1292     etaRange[1] = -999.;
1293
1294     return kFALSE; //FIXME: Shall we add something to invalidate this event in SetEvent to dont count it in AliAnalysisTaskMuMu::FillCounters?
1295   }
1296   else
1297   {
1298     etaRange[0] = etaMin;
1299     etaRange[1] = etaMax;
1300
1301     return kTRUE;
1302   }
1303 }
1304
1305 //_____________________________________________________________________________
1306 Double_t AliAnalysisMuMuNch::GetSPDCorrection(Double_t zvert, Double_t eta) const
1307 {
1308   if ( !fSPDOneOverAccxEff )
1309   {
1310     return 1.; // Like this we will simply compute the raw tracklets or if we are correctiong by the mean tracklets the eta bin will be counted as valid
1311   }
1312   else
1313   {
1314     Int_t bin =fSPDOneOverAccxEff->FindBin(zvert,eta);
1315     return fSPDOneOverAccxEff->GetBinContent(bin);// Like this we will compute Nch from the SPD AccxEFf
1316   }
1317 }
1318
1319 //_____________________________________________________________________________
1320 Double_t AliAnalysisMuMuNch::GetTrackletsMeanCorrection(Double_t zvert, Int_t nTracklets, Bool_t corrToCompare) const
1321 {
1322   if ( !fSPDMeanTracklets )
1323   {
1324     AliFatal("ERROR: No tracklets mean correction");
1325     return 0;
1326   }
1327   else
1328   {
1329     TProfile* h(0x0);
1330     if ( corrToCompare && fSPDMeanTrackletsCorrToCompare ) h = static_cast<TProfile*>(fSPDMeanTrackletsCorrToCompare);
1331     else h = static_cast<TProfile*>(fSPDMeanTracklets);
1332
1333     Int_t bin = h->FindBin(zvert);
1334     Double_t mNTtklsZ = h->GetBinContent(bin);
1335
1336     Double_t mNTtklsZref(0.);
1337     if ( fSPDMeanTrackletsCorrToCompare )
1338     {
1339       if ( corrToCompare ) mNTtklsZref = h->GetMaximum();
1340       else
1341       {
1342         if ( fMeanTrRef > 0. ) mNTtklsZref = fMeanTrRef;
1343         else mNTtklsZref = h->GetMaximum();
1344       }
1345     }
1346     else
1347     {
1348       if ( fMeanTrRef > 0. ) mNTtklsZref = fMeanTrRef;
1349       else mNTtklsZref = h->GetMaximum();
1350     }
1351
1352     Double_t deltaN;
1353     if ( mNTtklsZ == 0. ) deltaN = -1000.; // If the zvertex is out of the correction range the correction is < -999 (non valid)
1354     else
1355     {
1356       Double_t coef(1.);
1357       if ( mNTtklsZref < mNTtklsZ ) coef = -1.;
1358
1359       Double_t meanPoiss = nTracklets*(mNTtklsZref - mNTtklsZ)/mNTtklsZ;
1360       deltaN = coef*frand->PoissonD(TMath::Abs(meanPoiss));
1361     }
1362     return deltaN;
1363   }
1364 }
1365
1366
1367 //_____________________________________________________________________________
1368 AliAODTracklets* AliAnalysisMuMuNch::GetTracklets(AliVEvent* event)
1369 {
1370   /// Return AliAODTracklets if the event is an AOD event. If event is ESD event creates AliAODTracklets object from ESD AliMultiplicity and returns it
1371
1372   AliAODTracklets* tracklets(0x0);
1373
1374   if ( event->IsA() == AliAODEvent::Class() )
1375   {
1376     tracklets = static_cast<const AliAODEvent*>(Event())->GetTracklets(); // Tracklets object
1377   }
1378   else if ( event->IsA() == AliESDEvent::Class() )
1379   {
1380     const AliMultiplicity* mult = static_cast<const AliESDEvent*>(Event())->GetMultiplicity();
1381
1382     if (mult)
1383     {
1384       tracklets = new AliAODTracklets();
1385
1386       if ( mult->GetNumberOfTracklets() > 0 )
1387       {
1388         tracklets->CreateContainer(mult->GetNumberOfTracklets());
1389         for (Int_t n = 0 ; n < mult->GetNumberOfTracklets() ; n++)
1390         {
1391           tracklets->SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1392         }
1393       }
1394     }
1395   }
1396   else AliFatal("Unrecognized Event Type");
1397
1398   return tracklets;
1399
1400 }
1401
1402 //==============================These methods are useless here, anyway they should go in the EventCutterClass=================
1403 //_____________________________________________________________________________
1404 Bool_t AliAnalysisMuMuNch::HasAtLeastNTrackletsInEtaRange(const AliVEvent& event, Int_t n, Double_t& etaMin, Double_t& etaMax) const
1405 {
1406   if ( event.IsA() != AliAODEvent::Class() )
1407   {
1408     return kFALSE;
1409   }
1410
1411   AliAODTracklets* tracklets = static_cast<const AliAODEvent&>(event).GetTracklets();
1412
1413   if (!tracklets)
1414   {
1415     return kFALSE;
1416   }
1417
1418   Int_t nTrackletsInRange(0);
1419
1420   Int_t nTracklets = tracklets->GetNumberOfTracklets();
1421
1422   for (Int_t i = 0 ; i < nTracklets && nTrackletsInRange < n; i++)
1423   {
1424     Double_t eta = -TMath::Log(TMath::Tan(tracklets->GetTheta(i)/2.0));
1425
1426     if ( eta > etaMin && eta < etaMax )
1427     {
1428       ++nTrackletsInRange;
1429     }
1430   }
1431
1432   return (nTrackletsInRange>=n);
1433 }
1434
1435 //_____________________________________________________________________________
1436 void AliAnalysisMuMuNch::NameOfHasAtLeastNTrackletsInEtaRange(TString& name, Int_t n, Double_t& etaMin, Double_t& etaMax) const
1437 {
1438   if ( TMath::AreEqualAbs(TMath::Abs(etaMin),TMath::Abs(etaMax),1E-9 ) )
1439   {
1440     name.Form("ATLEAST%dTRKLINABSETALT%3.1f",n,TMath::Abs(etaMin));
1441   }
1442   else
1443   {
1444     name.Form("ATLEAST%dTRKLINETA%3.1f-%3.1f",n,etaMin,etaMax);
1445   }
1446 }
1447 //=======================================================================================
1448
1449 //_____________________________________________________________________________
1450 Bool_t AliAnalysisMuMuNch::IsMCtrackFromGenerator(Int_t indexMC) const
1451 {
1452   ///Checks if the MC particle corresponding to a track has been generated with a given generator.
1453
1454   if ( !HasMC() )
1455   {
1456     AliWarning("There is no MC information in the event");
1457     return kFALSE;
1458   }
1459
1460   if ( indexMC >= MCEvent()->GetNumberOfTracks() )
1461   {
1462     AliWarning("This MC track does not exist. Index out of range");
1463 //    return kFALSE;
1464   }
1465
1466 //______________"MANUAL" METHOD___________//
1467   //It does not work since NProduced gives just the stable dpmjet produced particles. See AliMCEvent::GetGenerator
1468 //  Int_t labelTemp(labelMC);
1469 //  while ( labelTemp != -1 )
1470 //  {
1471 //    AODpartMCTemp = static_cast<AliAODMCParticle*>(MCEvent()->GetTrack(labelMC));
1472 //    labelTemp = AODpartMCTemp->GetMother();
1473 //    if ( labelTemp != -1 ) labelMC = labelTemp;
1474 //    std::cout << labelTemp << " ===>" << labelMC <<std::endl;
1475 //  }
1476 //
1477 //  AliAODMCHeader* mcHeader = static_cast<AliAODMCHeader*>(Event()->FindListObject(AliAODMCHeader::StdBranchName()));
1478 //  if(mcHeader)
1479 //  {
1480 //    TList* lheaders = mcHeader->GetCocktailHeaders();
1481 //    lheaders->Print();
1482 //
1483 //    AliGenEventHeader* mcGenH(0x0);
1484 //    TIter next(lheaders); // Get the iterator on the list of cocktail headers
1485 //    Int_t nProduced(0);
1486 //    while ( labelMC >= nProduced ) // Loop over the cocktail headers
1487 //    {
1488 //      mcGenH = static_cast<AliGenEventHeader*>(next());
1489 //      if ( !mcGenH )
1490 //      {
1491 //        std::cout << "Generator header not found" << std::endl;
1492 //        return kFALSE;
1493 //      }
1494 //      else nProduced += mcGenH->NProduced();
1495 //    }
1496 //  }
1497 //  TString genName(mcGenH->GetName());
1498 //  return genName.Contains("dpmjet_0");
1499  //_________________________________________
1500
1501   TString genName;
1502   if ( !MCEvent()->GetCocktailGenerator(indexMC,genName) )
1503   {
1504     AliWarning("No cocktail generator found for this event");
1505     return kFALSE;
1506   }
1507 //  std::cout << genName.Data() << std::endl;
1508   return genName.Contains(fGeneratorHeaderClass->Data());
1509 }
1510
1511 //==============================This method is not used=====================
1512 //_____________________________________________________________________________
1513 Double_t AliAnalysisMuMuNch::NumberOfTrackletsInEtaRange(const AliVEvent& event, Double_t& etaMin, Double_t& etaMax, Bool_t /*corrected*/) const
1514 {
1515   if ( event.IsA() != AliAODEvent::Class() )
1516   {
1517     AliError("Not working for ESDs...");
1518     return 0;
1519   }
1520
1521   AliAODTracklets* tracklets = static_cast<const AliAODEvent&>(event).GetTracklets();
1522
1523   if (!tracklets)
1524   {
1525     return 0;
1526   }
1527
1528   Double_t nTrackletsInRange(0);
1529
1530   Int_t nTracklets = tracklets->GetNumberOfTracklets();
1531
1532   for (Int_t i = 0 ; i < nTracklets ; i++)
1533   {
1534     Double_t thetaTracklet = tracklets->GetTheta(i);
1535     Double_t etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
1536
1537     if ( etaTracklet > etaMin && etaTracklet < etaMax )
1538     {
1539       nTrackletsInRange += 1.0;
1540     }
1541   }
1542
1543   return nTrackletsInRange;
1544 }
1545 //==============================
1546
1547
1548 //_____________________________________________________________________________
1549 void AliAnalysisMuMuNch::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
1550 {
1551   /// Set the event, compute multiplicities and add them as a list to the event
1552
1553 //  if ( event->IsA() != AliAODEvent::Class() )
1554 //  {
1555 //    AliError("Only working for AODs for the moment.");
1556 //    return;
1557 //  }
1558
1559   AliAnalysisMuMuBase::SetEvent(event,mcEvent); // To have Event() and MCEvent() method working
1560
1561   TList* nchList = static_cast<TList*>(event->FindListObject("NCH")); // Define the list with the NCH info for each event
1562   if (!nchList)
1563   {
1564     nchList = new TList;
1565     nchList->SetOwner(kTRUE);
1566     nchList->SetName("NCH");
1567     event->AddObject(nchList);
1568   }
1569
1570   nchList->Clear(); // We clear the NCH list for this new event
1571
1572   Double_t SPDZv(0.);
1573   Bool_t vertexSPDFound(kFALSE);
1574 //  AliAODTracklets* tracklets(0x0);
1575
1576   AliVVertex* vertexSPD = AliAnalysisMuonUtility::GetVertexSPD(event);
1577   if ( vertexSPD )
1578   {
1579     vertexSPDFound = kTRUE;
1580     SPDZv = vertexSPD->GetZ();
1581     if ( SPDZv == 0. ) SPDZv = -40.; // If SPDZv = 0. means that there is no reconstructed vertex. Setting the vertex to -40 allow us to get an invalid correction
1582   }
1583
1584   AliAODTracklets* tracklets = GetTracklets(event);
1585
1586 //  if ( event->IsA() == AliAODEvent::Class() )
1587 //  {
1588 //    const AliAODVertex* vertexSPD = static_cast<const AliAODEvent*>(Event())->GetPrimaryVertexSPD(); // SPD vertex object
1589 //    if ( vertexSPD )
1590 //    {
1591 //      vertexSPDFound = kTRUE;
1592 //      SPDZv = vertexSPD->GetZ();
1593 //      if ( SPDZv == 0. ) SPDZv = -40.; // If SPDZv = 0. means that there is no reconstructed vertex. Setting the vertex to -40 allow us to get an invalid correction
1594 //    }
1595 ////    else return;
1596 //
1597 //    tracklets = static_cast<const AliAODEvent*>(Event())->GetTracklets(); // Tracklets object
1598 //  }
1599 //
1600 //  else if ( event->IsA() == AliESDEvent::Class() )
1601 //  {
1602 //    const AliESDVertex* vertexSPD = static_cast<const AliESDEvent*>(Event())->GetPrimaryVertexSPD(); // SPD vertex object
1603 //    if ( vertexSPD )
1604 //    {
1605 //      vertexSPDFound = kTRUE;
1606 //      SPDZv = vertexSPD->GetZ();
1607 //    }
1608 ////    else return;
1609 //
1610 //    tracklets = new AliAODTracklets();
1611 //    const AliMultiplicity* mult = static_cast<const AliESDEvent*>(Event())->GetMultiplicity();
1612 //
1613 //    if (mult)
1614 //    {
1615 //      if (mult->GetNumberOfTracklets()>0)
1616 //      {
1617 //        tracklets->CreateContainer(mult->GetNumberOfTracklets());
1618 //        for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++)
1619 //        {
1620 ////          if(fMChandler)
1621 ////          {
1622 ////            fMChandler->SelectParticle(mult->GetLabel(n, 0));
1623 ////            fMChandler->SelectParticle(mult->GetLabel(n, 1));
1624 ////          }
1625 //          tracklets->SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1626 //        }
1627 //      }
1628 //    }
1629 //  }
1630 //  else AliFatal("Unrecognized Event Type");
1631
1632   TH1* hSPDcorrectionVsEta(0x0); // Pointers for the individual event histos
1633   TH1* hNchVsEta(0x0);
1634   TH1* hNTrackletVsEta(0x0);
1635   TH1* hNTrackletVsPhi(0x0);
1636
1637   // Variables we will use in the multiplicity computation
1638   Int_t binMin,binMax;
1639   Int_t nTracklets(0);
1640   Double_t thetaTracklet(0.),phiTracklet(0.),etaTracklet(0.);
1641   Double_t nch(0.0);
1642   Int_t nTrackletsInRange(0);
1643   Int_t nTrackletsInRangeSPDOk(0);
1644   Int_t nBins(0);
1645
1646
1647   Bool_t isEtaRangeValid = GetEtaRangeSPD(SPDZv,fetaRange);
1648
1649   //____Data (or reco) multiplicity computation ___
1650   if ( !fResolution ) // When computing resolutions we dont do anything else
1651   {
1652     //_______Create once the histos with the individual event "properties" (cleared at the beginning of each new event)_______
1653     if ( !Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta") )
1654     {
1655       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","SPDcorrectionVsEta","SPD correction-like vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1656
1657       CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsEta","Nch vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1658
1659       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsEta","Ntracklet vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1660
1661       CreateEventHistos(kHistoForData,"AliAnalysisMuMuNch","NBkgTrackletsVSEta","NBkg Tracklets vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1662
1663       Double_t phimin = 0.; //Phi range
1664       Double_t phimax = 2*TMath::Pi();
1665       Int_t nphibins = GetNbins(phimin,phimax,0.05);
1666
1667       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsPhi","Ntracklet vs #phi;#phi",nphibins,phimin,phimax);
1668
1669       CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsPhi","Nch vs #phi;#phi",nphibins,phimin,phimax);
1670
1671       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","test","test",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1672
1673       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD1LR->Clone());
1674       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD1LL->Clone());
1675       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD2LR->Clone());
1676       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD2LL->Clone());
1677     }
1678     //_________
1679
1680     hSPDcorrectionVsEta = Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta"); // Set the individual event histos
1681     hNTrackletVsEta = Histo("AliAnalysisMuMuNch","NTrackletVsEta");
1682     hNTrackletVsPhi = Histo("AliAnalysisMuMuNch","NTrackletVsPhi");
1683
1684     hSPDcorrectionVsEta->Reset(); // Reset of the individual event histos
1685     hNTrackletVsEta->Reset();
1686     hNTrackletVsPhi->Reset();
1687
1688     //____Data (or Reco if simu) tracklets computation___
1689     if ( tracklets && vertexSPDFound && isEtaRangeValid )
1690     {
1691
1692       Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,fetaRange[1]);
1693       Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,fetaRange[0]);
1694
1695       nTracklets = tracklets->GetNumberOfTracklets();
1696
1697       Double_t SPDr; // SPD efficiency (If we use the 'mean correction' or no correction it will be always 1)
1698       for (Int_t i = 0 ; i < nTracklets ; i++)
1699       {
1700         thetaTracklet = tracklets->GetTheta(i);
1701         etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
1702         if ( etaTracklet < fetaRange[0] || etaTracklet > fetaRange[1] ) continue; // Avoid tracklets out of the eta SPD acceptance or out of the eta cut
1703
1704         SPDr = GetSPDCorrection(SPDZv,etaTracklet); // Get SPD AccxEff for the corresponding [zvtx,eta] bin (If we use the 'mean correction' or no correction it will be always 1)
1705
1706         Int_t bin = fEtaAxis->FindBin(etaTracklet);
1707
1708         nTrackletsInRange++;
1709         if ( SPDr!=0. && SPDr <= 2.5) // Threshold to reduce border effects
1710         {
1711           hSPDcorrectionVsEta->SetBinContent(bin,SPDr);
1712           hNTrackletVsEta->Fill(etaTracklet);
1713           hNTrackletVsPhi->Fill(tracklets->GetPhi(i));
1714
1715           nTrackletsInRangeSPDOk++;
1716         }
1717         else // If the correction is above the threshold we store a -1 in the correction to skip this eta bin int the fill method (can just happen using SPD AccxEff correction)
1718         {
1719           hSPDcorrectionVsEta->SetBinContent(bin,-1.0);
1720         }
1721       }
1722
1723       nchList->Add(new TParameter<Double_t>("Ntr",nTrackletsInRange));
1724       nchList->Add(new TParameter<Double_t>("NtrSPDOk",nTrackletsInRangeSPDOk)); // if there is not SPD AccxEff correction this will be equal to Ntr
1725
1726       //___ Fill the out-of-eta-range bins with -1.0
1727
1728       binMin = fEtaAxis->FindBin(fetaRange[0]);
1729       binMax = fEtaAxis->FindBin(fetaRange[1]);
1730
1731       for ( Int_t i = 1; i < binMin; ++i )
1732       {
1733         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
1734       }
1735       for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
1736       {
1737         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
1738       }
1739
1740       if ( fSPDOneOverAccxEff || fSPDMeanTracklets ) // Here we correct the raw tracklets
1741       {
1742         if ( fSPDOneOverAccxEff ) // In this case the correction is the SPD accxEff
1743         {
1744
1745           //----- Mean dNchdEta computation to set it into the event
1746           for (Int_t j = binMin/*1*/ ; j <= binMax/*fEtaAxis->GetNbins()*/ ; j++) // Loop over eta bins
1747           {
1748             Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
1749
1750             Double_t eta = fEtaAxis->GetBinCenter(j);
1751
1752             if ( correction < 0 ) continue; // If the correction is < 0 we skip that eta bin
1753             else if ( correction == 0.0 ) // If the correction is 0 we have no tracklets in that eta bin
1754             {
1755               Double_t spdCorr = GetSPDCorrection(SPDZv,eta);
1756               if ( spdCorr == 0. || spdCorr > 2.5) continue; // If the correction in the eta bin is not within the threshold we do not count the "0"(that eta bin will not count for nBins)
1757             }
1758
1759             nch += hNTrackletVsEta->GetBinContent(j) * correction; // Number of charged particles (tracklets*(1/SPDAccxEff))
1760
1761             ++nBins; // We sum up the number of bins entering in the computation
1762           }
1763
1764           nchList->Add(new TParameter<Double_t>("Nch",nch));
1765
1766           Double_t meandNchdEta(0.);// We compute the mean dNch/dEta in the event
1767           if ( nBins >  0 )
1768           {
1769             meandNchdEta = nch / (nBins*fEtaAxis->GetBinWidth(5));
1770           }
1771
1772           nchList->Add(new TParameter<Double_t>("MeandNchdEta",meandNchdEta)); // We add the mean dNch/dEta to the event. It will serve us as a multiplicity estimator.
1773         }
1774 //        else if( fSPDMeanTracklets ) // In this case the correction is the 'mean correction'
1775         if( fSPDMeanTracklets ) // In this case the correction is the 'mean correction'
1776         {
1777           SPDr = GetTrackletsMeanCorrection(SPDZv,nTrackletsInRange); // Get 'mean correction' for the zvtx
1778
1779           if ( SPDr < -999.) nch = -1;
1780           else nch = nTrackletsInRange + SPDr;
1781
1782           nchList->Add(new TParameter<Double_t>("NtrCorr",nch));// We add the corrected number of tracklets to the event. It will serve us as a multiplicity estimator.
1783         }
1784 //        else AliFatal("Unknown tracklets correction class");
1785       }
1786
1787     }
1788     else // If there is no tracklets or vertex object or the eta range is not valid everything is invalidated
1789     {
1790       for ( Int_t i = 1 ; i <= hSPDcorrectionVsEta->GetNbinsX() ; ++i )
1791       {
1792         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
1793       }
1794
1795       nchList->Add(new TParameter<Double_t>("Ntr",-1)); //We set an invalid Ntracklets
1796
1797       if ( fSPDOneOverAccxEff ) // In this case the correction is the SPD accxEff
1798       {
1799         nchList->Add(new TParameter<Double_t>("MeandNchdEta",-1)); //We set an invalid dNch/deta
1800       }
1801       else if( fSPDMeanTracklets ) // In this case the correction is the 'mean correction'
1802       {
1803         nchList->Add(new TParameter<Double_t>("NtrCorr",-1)); //We set an invalid Ntracklets corrected
1804       }
1805     }
1806     //_______
1807
1808     //____ We add the tracklets histos to the event
1809     nchList->Add(hSPDcorrectionVsEta->Clone());
1810     nchList->Add(hNTrackletVsEta->Clone());
1811     nchList->Add(hNTrackletVsPhi->Clone());
1812     //_______
1813   }
1814   //_______
1815
1816   //____Input MC multiplicity computation ___
1817   if ( HasMC() )
1818   {
1819     Double_t etaRange[2];
1820     if ( !fResolution ) //When computing resolutions we dont do anything else
1821     {
1822       Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent());
1823       GetEtaRangeSPD(MCZv,etaRange);
1824
1825       hNchVsEta = MCHisto("AliAnalysisMuMuNch","NchVsEta");
1826       hSPDcorrectionVsEta = MCHisto("AliAnalysisMuMuNch","SPDcorrectionVsEta");
1827       TH1* hNchVsPhi = MCHisto("AliAnalysisMuMuNch","NchVsPhi");
1828       hNTrackletVsEta = MCHisto("AliAnalysisMuMuNch","NTrackletVsEta");
1829       hNTrackletVsPhi = MCHisto("AliAnalysisMuMuNch","NTrackletVsPhi");
1830
1831       hNTrackletVsEta->Reset();
1832       hNchVsEta->Reset();
1833       hNchVsPhi->Reset();
1834       hSPDcorrectionVsEta->Reset();
1835       hNTrackletVsPhi->Reset();
1836
1837       //___ Fill the out-of-eta-range bins with -1.0
1838       binMin = fEtaAxis->FindBin(etaRange[0]);
1839       binMax = fEtaAxis->FindBin(etaRange[1]);
1840
1841       for ( Int_t i = 1; i < binMin; ++i )
1842       {
1843         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
1844       }
1845       for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
1846       {
1847         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
1848       }
1849       for ( Int_t i = binMin; i <= binMax; ++i ) // Fill the bins inside the eta range with +1
1850       {
1851         hSPDcorrectionVsEta->SetBinContent(i,1.0); // FIXME: Do we really want all the tracklets to be counted for MC purposes? or we want to count only those in active regions of the SPD as for Data computations?
1852       }
1853       //___
1854
1855       Int_t nMCTracks = MCEvent()->GetNumberOfTracks(); // number of MC tracks
1856
1857       for ( Int_t i = 0; i < nMCTracks ; ++i ) //Loop over generated tracks
1858       {
1859 //        AliAODMCParticle* AODpart = static_cast<AliAODMCParticle*>(mcEvent->GetTrack(i));
1860         AliMCParticle* MCpart = static_cast<AliMCParticle*>(mcEvent->GetTrack(i));
1861
1862         Bool_t isPP(kFALSE);
1863         Int_t nGentorDex(-2); //Generator Index: -1 as default;
1864         //                  0 is for DPMJET;
1865         //                  1, 2 ... are for other generators
1866
1867         if ( static_cast<const AliVEvent*>(event)->IsA() == AliAODEvent::Class() )
1868         {
1869           AliAODMCParticle* AODpart = static_cast<AliAODMCParticle*>(mcEvent->GetTrack(i));
1870           isPP = AODpart->IsPhysicalPrimary();
1871           nGentorDex = AODpart->GetGeneratorIndex();
1872         }
1873         else if ( static_cast<const AliVEvent*>(event)->IsA() == AliESDEvent::Class() )
1874         {
1875           isPP = MCEvent()->IsPhysicalPrimary(i);
1876           nGentorDex = MCpart->GetGeneratorIndex();
1877         }
1878
1879
1880 //        if ( AODpart->IsPhysicalPrimary() ) // We take only particles produced in the collision (Particles produced in the collision including products of strong and electromagnetic decay and excluding feed-down from weak decays of strange particles)
1881         if ( isPP && nGentorDex==0 )
1882         {
1883 //          if ( !IsMCtrackFromGenerator(i) ) continue; // Select only the particles generated by the desired generator
1884
1885           if ( MCpart->Charge()!=0 ) // We take only charged particles
1886           {
1887             hNchVsEta->Fill(MCpart->Eta());
1888             hNchVsPhi->Fill(MCpart->Phi());
1889
1890           }
1891         }
1892       }
1893
1894       nchList->Add(hNchVsEta->Clone("MCNchVsEta"));
1895       nchList->Add(hNchVsPhi->Clone("MCNchVsPhi"));
1896       nchList->Add(hSPDcorrectionVsEta->Clone("MCSPDcorrectionVsEta"));
1897     }
1898               //__Bkg tracklets and Resolution estimation __
1899     if ( tracklets ) // We can compute the Bkg tracklets and resolution only if we have the tracklets object in the event
1900     {
1901       TH1* hNBkgTrackletsVSEta(0x0); // Pointer for the Bkg histo
1902       if ( !fResolution )
1903       {
1904         hNBkgTrackletsVSEta = Histo("AliAnalysisMuMuNch","NBkgTrackletsVSEta");
1905
1906         hNBkgTrackletsVSEta->Reset();
1907       }
1908
1909       nTracklets = tracklets->GetNumberOfTracklets();
1910
1911       for (Int_t i = 0 ; i < nTracklets ; i++) // Loop over tracklets to check if they come or not from the same MC particle
1912       {
1913         Int_t label1 = tracklets->GetLabel(i,0);
1914         Int_t label2 = tracklets->GetLabel(i,1);
1915
1916 //        if ( !IsMCtrackFromGenerator(label1) || !IsMCtrackFromGenerator(label2) ) continue; // Select only the particles generated by the desired generator
1917
1918         thetaTracklet = tracklets->GetTheta(i);
1919         etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
1920         phiTracklet = tracklets->GetPhi(i);
1921
1922         if ( !fResolution )
1923         {
1924           hNTrackletVsEta->Fill(etaTracklet);
1925           hNTrackletVsPhi->Fill(phiTracklet);
1926         }
1927
1928         if ( label1 != label2 && !fResolution ) hNBkgTrackletsVSEta->Fill(etaTracklet); // Tracklets not comming from the same MC particle are Bkg
1929         else if ( !fSPDOneOverAccxEff && !fSPDMeanTracklets && fResolution ) // Compute the resolutions with the tracklets comming from the same MC particle
1930         {
1931 //          AliAODMCParticle* AODpartMC = static_cast<AliAODMCParticle*>(MCEvent()->GetTrack(label1));
1932           AliMCParticle* partMC = static_cast<AliMCParticle*>(MCEvent()->GetTrack(label1));
1933           Double_t etaTrackletMC = partMC->Eta();
1934           Double_t phiTrackletMC = partMC->Phi();
1935
1936           // Resolution variables
1937           nchList->Add(new TParameter<Double_t>(Form("EtaReco%d",label1),etaTracklet));
1938           nchList->Add(new TParameter<Double_t>(Form("EtaMC%d",label1),etaTrackletMC));
1939
1940           nchList->Add(new TParameter<Double_t>(Form("PhiReco%d",label1),phiTracklet));
1941           nchList->Add(new TParameter<Double_t>(Form("PhiMC%d",label1),phiTrackletMC));
1942
1943         }
1944       }
1945
1946       if (!fResolution )
1947       {
1948         nchList->Add(hNTrackletVsEta->Clone("MCNTrackletVsEta"));
1949         nchList->Add(hNTrackletVsPhi->Clone("MCNTrackletVsPhi"));
1950         nchList->Add(hNBkgTrackletsVSEta->Clone());
1951       }
1952
1953     }
1954
1955   }
1956   //_______
1957
1958 }
1959
1960 //_____________________________________________________________________________
1961 void AliAnalysisMuMuNch::SetRun(const AliInputEventHandler* eventHandler)
1962 {
1963   // For each new run this method sets the corresponding SPD correction from SPDCorrection list using the SPDCorrection map
1964
1965   if ( !fSPDCorrectionList ) return;
1966
1967 //  AliAODEvent* event = static_cast<AliAODEvent*>(eventHandler->GetEvent());
1968    AliVEvent* event = static_cast<AliVEvent*>(eventHandler->GetEvent());
1969
1970   TString run(Form("%d",event->GetRunNumber())); // Get the run number
1971   TObjString*  SPDCorrectionKey = static_cast<TObjString*>(fSPDCorrectionMap->GetValue(run)); // Get the corresponding SPDCorrection map key for the run
1972
1973   Int_t SPDCorrectionIndex(0);
1974   if ( SPDCorrectionKey ) SPDCorrectionIndex = SPDCorrectionKey->String().Atoi(); // Converts the key into Int, if found
1975   else AliFatal(Form("No SPD correction found for run %d",run.Atoi()));
1976
1977   TH1* hSPDCorrection = static_cast<TH1*>(fSPDCorrectionList->At(SPDCorrectionIndex)); // Gets the SPD correction at key position from the list
1978
1979   if ( static_cast<TH1*>(hSPDCorrection)->IsA() == TH2::Class() ) fSPDOneOverAccxEff = static_cast<TH2F*>(hSPDCorrection); // Sets the corresponding correction data member
1980   else if ( static_cast<TH1*>(hSPDCorrection)->IsA() == TProfile::Class() ) fSPDMeanTracklets = static_cast<TProfile*>(hSPDCorrection);
1981   else AliFatal("Unrecognized correction class");
1982
1983   AliInfo(Form("Using correction %s for run %s",hSPDCorrection->GetName(),run.Data()));
1984
1985 }
1986
1987 //_____________________________________________________________________________
1988 void AliAnalysisMuMuNch::Terminate(Option_t *)
1989 {
1990   /// Called once at the end of the query
1991   if ( !HistogramCollection() ) return;
1992
1993   if ( HistogramCollection()->FindObject(Form("/%s/AliAnalysisMuMuNch/NTrackletVsEta",MCInputPrefix())) )
1994   {
1995     HistogramCollection()->Remove(Form("/%s/AliAnalysisMuMuNch/NTrackletVsEta",MCInputPrefix()));
1996     HistogramCollection()->Remove(Form("/%s/AliAnalysisMuMuNch/NTrackletVsPhi",MCInputPrefix()));
1997     HistogramCollection()->Remove(Form("/%s/AliAnalysisMuMuNch/NchVsEta",MCInputPrefix()));
1998     HistogramCollection()->Remove(Form("/%s/AliAnalysisMuMuNch/NchVsPhi",MCInputPrefix()));
1999     HistogramCollection()->Remove(Form("/%s/AliAnalysisMuMuNch/SPDcorrectionVsEta",MCInputPrefix()));
2000     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NBkgTrackletsVSEta");
2001   }
2002
2003   if ( HistogramCollection()->FindObject("/AliAnalysisMuMuNch/NTrackletVsEta") )
2004   {
2005     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsEta");
2006     HistogramCollection()->Remove("/AliAnalysisMuMuNch/test");
2007     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsPhi");
2008     HistogramCollection()->Remove("/AliAnalysisMuMuNch/SPDcorrectionVsEta");
2009   }
2010   //____ Compute dNchdEta histo
2011   TObjArray* idArr =  HistogramCollection()->SortAllIdentifiers();
2012
2013   TIter next(idArr);
2014   TObjString* id;
2015
2016   while ( (id = static_cast<TObjString*>(next())) )
2017   {
2018     TProfile* p = static_cast<TProfile*>(HistogramCollection()->FindObject(Form("%s%s",id->GetName(),"MeanNchVsEta")));
2019
2020     if ( !p ) continue;
2021
2022     TH1* h = new TH1F("MeandNchdEta","Event averaged dN_{ch}/d#eta ;#eta;<dN_{ch}/d#eta>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
2023
2024     if ( p->GetNbinsX() != h->GetNbinsX() || p->GetXaxis()->GetXmin() != h->GetXaxis()->GetXmin() || p->GetXaxis()->GetXmax() != h->GetXaxis()->GetXmax() )
2025     {
2026       AliError("ERROR: Cannot compute MeandNchdEta since the binning doesn't match with MeanNchVsEta histo");
2027       continue;
2028     }
2029
2030     for ( Int_t i = 1 ; i < h->GetNbinsX() ; i++ )
2031     {
2032       h->SetBinContent(i,p->GetBinContent(i)/p->GetBinWidth(i));
2033       h->SetBinError(i,p->GetBinError(i)/p->GetBinWidth(i));
2034       h->SetEntries(p->GetEntries());
2035     }
2036
2037     HistogramCollection()->Adopt(Form("%s",id->GetName()),h);
2038   }
2039
2040   delete idArr;
2041   //___
2042 }