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 "AliAODTracklets.h"
  18 #include "AliAnalysisMuonUtility.h"
  19 #include "AliLog.h"
  20 #include "AliAODEvent.h"
  21 #include "AliESDUtils.h"
  22 #include "TMath.h"
  23 #include "AliAnalysisMuMuCutRegistry.h"
  24 #include "AliAnalysisMuMuCutElement.h"
  25 #include "Riostream.h"
  26 #include "TParameter.h"
  27 #include <set>
  28 #include <utility>
  29 #include "AliMergeableCollection.h"
  30 #include "AliMCEvent.h"
  31 #include "AliAODMCParticle.h"
  32 #include "TAxis.h"
  33 #include "TCanvas.h"
  34 #include "TH2F.h"
  35 #include "TF1.h"
  36 #include "TProfile.h"
  37 #include "TObjString.h"
  38
  39
  40 namespace {
  41
  42   Double_t SPDgeomR(Double_t* x,Double_t* par) // Eta position of the SPD right edge eta as "seen" from a z vertex position
  43   {
  44     // par[0] = radius of SPD layer
  45
  46     Double_t d = x[0];
  47     Double_t z(0.);
  48     Double_t theta(0.);
  49
  50     if( d < 14.09999 )
  51     {
  52       z = 14.1 - d;
  53       theta = TMath::ATan(par[0]/z);
  54     }
  55
  56     else if( d > 14.09999 )
  57     {
  58       z = d - 14.1;
  59       theta = TMath::Pi() - TMath::ATan(par[0]/z);
  60     }
  61
  62     return -TMath::Log(TMath::Tan(theta/2.));
  63   }
  64
  65   Double_t SPDgeomL(Double_t* x,Double_t* par) // Eta position of the SPD left edge eta as "seen" from a z vertex position
  66   {
  67     // par[0] = radius of SPD layer
  68
  69     Double_t d = x[0];
  70     Double_t z(0.);
  71     Double_t theta(0.);
  72
  73     if( d > -14.09999 )
  74     {
  75       z = 14.1 + d;
  76       theta = TMath::Pi() - TMath::ATan(par[0]/z);
  77     }
  78
  79     if( d < -14.09999 )
  80     {
  81       z = -14.1 - d;
  82       theta = TMath::ATan(par[0]/z);
  83     }
  84
  85     return -TMath::Log(TMath::Tan(theta/2.));
  86   }
  87
  88 }
  89
  90 //_____________________________________________________________________________
  91 AliAnalysisMuMuNch::AliAnalysisMuMuNch(TH2* spdCorrection, Double_t etaMin, Double_t etaMax
  92                                       , Double_t zMin, Double_t zMax,Bool_t disableHistos, Bool_t computeResolution)
  93 : AliAnalysisMuMuBase(),
  94 fSPDCorrection(0x0),
  95 fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
  96 fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
  97 fCurrentEvent(0x0),
  98 fEtaMin(etaMin),
  99 fEtaMax(etaMax),
 100 fZMin(zMin),
 101 fZMax(zMax),
 102 fResolution(computeResolution)
 103 {
 104   if ( spdCorrection )
 105   {
 106     fSPDCorrection = static_cast<TH2F*>(spdCorrection->Clone());
 107     fSPDCorrection->SetDirectory(0);
 108   }
 109   DefineSPDAcceptance();
 110
 111   if ( disableHistos )
 112   {
 113     DisableHistograms("*");
 114   }
 115 }
 116
 117 //_____________________________________________________________________________
 118 AliAnalysisMuMuNch::~AliAnalysisMuMuNch()
 119 {
 120   delete fSPDCorrection;
 121   delete fEtaAxis;
 122   delete fZAxis;
 123   delete fSPD1LR;
 124   delete fSPD1LL;
 125   delete fSPD2LR;
 126   delete fSPD2LL;
 127 }
 128
 129 //_____________________________________________________________________________
 130 void AliAnalysisMuMuNch::DefineHistogramCollection(const char* eventSelection,
 131                                                    const char* triggerClassName,
 132                                                    const char* centrality)
 133 {
 134   // Define multiplicity histos
 135
 136   if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuNch") )
 137   {
 138     return;
 139   }
 140
 141   // dummy histogram to signal that we already defined all our histograms (see above)
 142   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"AliAnalysisMuMuNch","Dummy semaphore",1,0,1);
 143
 144   Double_t multMin = 0.;  //Tracklets multiplicity range
 145   Double_t multMax = 500.;
 146   Int_t nbinsMult = GetNbins(multMin,multMax,1.);
 147
 148   Double_t phimin = 0.; //Phi range
 149   Double_t phimax = 2*TMath::Pi();
 150   Int_t nphibins = GetNbins(phimin,phimax,0.05);
 151
 152   if ( !fSPDCorrection && fResolution ) // Resolution histos
 153   {
 154     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"EtaRes","#eta resolution;#eta_{Reco} - #eta_{MC};Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 155
 156     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiRes","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*2,-phimax,phimax);
 157
 158     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiResShifted","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*4,-phimax/4,phimax/4);
 159
 160     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());
 161
 162     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);
 163
 164     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());
 165
 166     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);
 167
 168     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());
 169
 170     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());
 171
 172     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Phi","Reco #phi distribution;#phi;Counts",nphibins,phimin,phimax);
 173
 174     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCPhi","MC #phi distribution;#phi;Counts",nphibins,phimin,phimax);
 175
 176     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","Reco #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 177
 178     CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCEta","MC #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 179
 180     return; // When computing resolutions we don't want to create the rest of histos
 181   }
 182
 183   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);
 184   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi");// Attach the SPD acc curves
 185
 186   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());
 187   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta");
 188
 189   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Tracklets","Number of tracklets distribution;N_{Tracklets};N_{events}",nbinsMult,multMin,multMax);
 190
 191   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());
 192   AttachSPDAcceptance(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NBkgTrackletsVsZVertexVsEta");
 193
 194   CreateEventHistos(kHistoForData | 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());
 195   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta");
 196
 197   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"NchVsZVertexVsPhi","Number of charged particles vs ZVertex vs #phi;ZVertex;#eta",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),nphibins,phimin,phimax);
 198
 199   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)
 200   AttachSPDAcceptance(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta");
 201
 202   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Nch","Number of charged particles distribution;N_{ch};N_{events}",nbinsMult,multMin,multMax);
 203
 204   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdEta","<dNchdEta> distribution;dN_{ch}/d#eta;N_{events}",nbinsMult,multMin,multMax);
 205
 206   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"TrackletsVsNch","Number of tracklets vs Number of charged particles;N_{ch};N_{tracklets}",nbinsMult,multMin,multMax,nbinsMult,multMin,multMax); //Response matrix
 207
 208   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"V0AMultVsNch","V0A multiplicity vs Number of charged particles;N_{ch};V0A Mult",nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 209
 210   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"V0CMultVsNch","V0C multiplicity vs Number of charged particles;N_{ch};V0C Mult",nbinsMult,multMin,multMax,nbinsMult,multMin,multMax);
 211
 212   // profile histograms
 213
 214   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta","Mean number of tracklets vs #eta;#eta;<N_{Tracklets}>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0);
 215
 216   CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsVsZVertex","Mean number of tracklets vs Z vertex;Z vertex;<N_{Tracklets}>",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 217
 218   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"MeanNchVsEta","Mean number of 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)
 219
 220    CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"MeanNchVsZVertex","Mean number of charged particles vs Z vertex;Z vertex;<N_{ch}>",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),0);
 221
 222 }
 223
 224 //_____________________________________________________________________________
 225 void AliAnalysisMuMuNch::DefineSPDAcceptance()
 226 {
 227   // Defines the functions ( eta = f(z) )of the edges (right/left) of the inner and outer SPD layers
 228   // R_in = 3.9 cm ; R_out = 7.6 cm
 229   fSPD1LR = new TF1("fSPD1LR",SPDgeomR,-40,40,1);
 230   fSPD1LR->SetParameter(0,3.9);
 231   fSPD1LL = new TF1("fSPD1LL",SPDgeomL,-40,40,1);
 232   fSPD1LL->SetParameter(0,3.9);
 233   fSPD2LR = new TF1("fSPD2LR",SPDgeomR,-40,40,1);
 234   fSPD2LR->SetParameter(0,7.6);
 235   fSPD2LL = new TF1("fSPD2LL",SPDgeomL,-40,40,1);
 236   fSPD2LL->SetParameter(0,7.6);
 237
 238 }
 239
 240 //_____________________________________________________________________________
 241 void AliAnalysisMuMuNch::AddHisto(const char* eventSelection,
 242                                   const char* triggerClassName,
 243                                   const char* centrality,
 244                                   const char* histoname,
 245                                   Double_t z,
 246                                   TH1* h,Bool_t isMC)
 247 {
 248   // Adds the content of a 1D histo to a 2D histo a the z position
 249
 250   Int_t zbin = fZAxis->FindBin(z);
 251   TH2F* h2;
 252   if (isMC) h2 = static_cast<TH2F*>(MCHisto(eventSelection,triggerClassName,centrality,histoname));
 253   else h2 = static_cast<TH2F*>(Histo(eventSelection,triggerClassName,centrality,histoname));
 254
 255   for ( Int_t i = 1; i <= h->GetXaxis()->GetNbins(); ++i )
 256   {
 257     Double_t content = h2->GetCellContent(zbin,i);
 258
 259     if ( h->GetBinContent(i) >  0 )
 260     {
 261       h2->SetCellContent(zbin,i,content + h->GetBinContent(i));
 262     }
 263   }
 264
 265   h2->SetEntries(h2->GetSumOfWeights());
 266 }
 267
 268
 269 //_____________________________________________________________________________
 270 void AliAnalysisMuMuNch::AttachSPDAcceptance(UInt_t dataType,
 271                                              const char* eventSelection,
 272                                              const char* triggerClassName,
 273                                              const char* centrality,const char* histoname)
 274 {
 275   if ( dataType & kHistoForData )
 276   {
 277     if( !Histo(eventSelection,triggerClassName,centrality,histoname) )
 278     {
 279       AliError(Form("ERROR: SPD Acceptance curves attach failed. Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
 280       return;
 281     }
 282
 283     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LR);
 284     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LL);
 285     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LR);
 286     Histo(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LL);
 287   }
 288   if ( dataType & kHistoForMCInput )
 289   {
 290     if( !MCHisto(eventSelection,triggerClassName,centrality,histoname) )
 291     {
 292       AliError(Form("ERROR: SPD Acceptance curves attach failed. MC Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
 293       return;
 294     }
 295
 296     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LR);
 297     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD1LL);
 298     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LR);
 299     MCHisto(eventSelection,triggerClassName,centrality,histoname)->GetListOfFunctions()->Add(fSPD2LL);
 300   }
 301
 302 }
 303
 304
 305 //_____________________________________________________________________________
 306 void AliAnalysisMuMuNch::FillHistosForEvent(const char* eventSelection,
 307                                             const char* triggerClassName,
 308                                             const char* centrality)
 309 {
 310   // Fills the data (or reco if simu) multiplicity histos
 311
 312   if ( IsHistogrammingDisabled() ) return;
 313
 314   if ( fResolution ) return; //When computing resolutions we skip this method
 315
 316   if ( !AliAnalysisMuonUtility::IsAODEvent(Event()) )
 317   {
 318     AliError("Don't know how to deal with ESDs...");
 319     return;
 320   }
 321
 322   if ( HasMC() && !fSPDCorrection ) // We have MC but no correction (SPD correction computation mode)so we skip the method
 323   {
 324     return;
 325   }
 326
 327   AliAODEvent* aod = static_cast<AliAODEvent*>(Event());
 328
 329   AliVVertex* vertex = aod->GetPrimaryVertexSPD();
 330
 331   TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
 332
 333   if (!nchList || nchList->IsEmpty() ) // Empty NCH means that there is no SPD vertex ( see SetEvent() ) when runing on data.
 334   {
 335     return;
 336   }
 337
 338   Double_t SPDZv;
 339
 340   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
 341   {
 342     SPDZv = -40.;
 343   }
 344
 345   else SPDZv = vertex->GetZ();
 346
 347   TH1* hSPDcorrectionVsEta = static_cast<TH1*>(nchList->FindObject("SPDcorrectionVsEta"));
 348   TH1* hNTrackletVsEta = static_cast<TH1*>(nchList->FindObject("NTrackletVsEta"));
 349   TH1* hNTrackletVsPhi = static_cast<TH1*>(nchList->FindObject("NTrackletVsPhi"));
 350
 351   TH1* hNchVsEta =static_cast<TH1*>(hNTrackletVsEta->Clone("NchVsEta"));
 352
 353   TProfile* hMeanTrackletsVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta"));
 354   TProfile* hMeanNchVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
 355
 356   TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(Histo(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
 357
 358   Int_t nBins(0);
 359
 360   Double_t nch(0.0);
 361   Double_t nTracklets(0.0);
 362
 363   for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) // Loop over eta bins
 364   {
 365     Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
 366
 367     Double_t eta = fEtaAxis->GetBinCenter(j);
 368
 369     if ( correction < 0 ) continue;
 370     else if ( correction == 0.0 ) // No tracklets found in this eta bin.
 371     {
 372       correction = GetSPDCorrection(SPDZv,eta);
 373
 374       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
 375     }
 376
 377     Double_t ntr = hNTrackletVsEta->GetBinContent(j); // Tracklets in eta bin
 378
 379     nch += ntr * correction; // Number of charged particles (corrected tracklets)
 380
 381     nTracklets += ntr;
 382
 383     hMeanTrackletsVsEta->Fill(eta,ntr); // Fill the number of tracklets of each eta bin in the profile
 384
 385     hMeanNchVsEta->Fill(eta,ntr*correction);
 386
 387     ++nBins; // We sum up the number of bins entering in the computation
 388
 389     // Fill the number of charged particles of each eta bin in the profile
 390     hNchVsEta->SetBinContent( j, hNchVsEta->GetBinContent(j) * correction );
 391
 392     hEventsVsZVertexVsEta->Fill(SPDZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
 393   }
 394
 395   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi",SPDZv,hNTrackletVsPhi);
 396   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta",SPDZv,hNTrackletVsEta);
 397
 398   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",SPDZv,hNchVsEta);
 399
 400   Histo(eventSelection,triggerClassName,centrality,"Tracklets")->Fill(nTracklets);
 401   Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsVsZVertex")->Fill(SPDZv,nTracklets);
 402   Histo(eventSelection,triggerClassName,centrality,"MeanNchVsZVertex")->Fill(SPDZv,nch);
 403
 404   Histo(eventSelection,triggerClassName,centrality,"TrackletsVsNch")->Fill(nTracklets,nch);
 405   Histo(eventSelection,triggerClassName,centrality,"Nch")->Fill(nch);
 406
 407   Double_t V0AMult = 0.;
 408   Double_t V0CMult = 0.;
 409
 410   AliVVZERO* vzero = aod->GetVZEROData();
 411   if (vzero)
 412   {
 413     Double_t multV0A = vzero->GetMTotV0A();
 414     V0AMult = AliESDUtils::GetCorrV0A(multV0A,SPDZv);
 415     Double_t multV0C = vzero->GetMTotV0C();
 416     V0CMult = AliESDUtils::GetCorrV0C(multV0C,SPDZv);
 417
 418     Histo(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nch);
 419     Histo(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nch);
 420
 421   }
 422
 423
 424   delete hNchVsEta;
 425
 426   // Mean dNch/dEta computation
 427   Double_t meandNchdEta(0.);
 428
 429   if ( nBins >  0 )
 430   {
 431     meandNchdEta = nch / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
 432   }
 433
 434   Histo(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
 435
 436
 437
 438 }
 439
 440 //_____________________________________________________________________________
 441 void AliAnalysisMuMuNch::FillHistosForMCEvent(const char* eventSelection,const char* triggerClassName,const char* centrality)
 442 {
 443   /// Fill input MC multiplicity histos
 444
 445   if ( IsHistogrammingDisabled() ) return;
 446
 447   TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
 448
 449   if (!nchList || nchList->IsEmpty())
 450   {
 451     return;
 452   }
 453
 454   Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent()); // Definition of MC generated z vertex
 455   AliVVertex* vertex = static_cast<AliAODEvent*>(Event())->GetPrimaryVertexSPD();
 456
 457   Double_t SPDZv(0.);
 458
 459   //____Resolution Histos___
 460   if ( !fSPDCorrection && fResolution )
 461   {
 462     Int_t nContributors(0);
 463     if ( vertex )
 464     {
 465       nContributors  = vertex->GetNContributors();
 466       SPDZv = vertex->GetZ();
 467     }
 468     MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsnC")->Fill(nContributors,SPDZv - MCZv);
 469     MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsMCz")->Fill(MCZv,SPDZv - MCZv);
 470
 471     Double_t EtaReco(0.),EtaMC(0.),PhiReco(0.),PhiMC(0.);
 472     Int_t i(-1),labelEtaReco(-1),labelEtaMC(-1),labelPhiReco(-1),labelPhiMC(-1);
 473
 474     while ( i < nchList->GetEntries() - 1 )
 475     {
 476       i++;
 477       while ( nchList->At(i)->IsA() != TObjString::Class() ) // In case there is a diferent object, just to skip it
 478       {
 479         i++;
 480       }
 481
 482       TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(nchList->At(i));
 483
 484
 485          //Eta Resolution
 486       if ( TString(p->GetName()).Contains("EtaReco") ) // We take the reco eta
 487       {
 488         sscanf(p->GetName(),"EtaReco%d",&labelEtaReco);
 489         EtaReco = p->GetVal();
 490         MCHisto(eventSelection,triggerClassName,centrality,"Eta")->Fill(EtaReco);
 491       }
 492       else if ( TString(p->GetName()).Contains("EtaMC") ) // We take the generated eta
 493       {
 494         sscanf(p->GetName(),"EtaMC%d",&labelEtaMC);
 495         EtaMC = p->GetVal();
 496         MCHisto(eventSelection,triggerClassName,centrality,"MCEta")->Fill(EtaMC);
 497       }
 498       if ( labelEtaReco > 0 && labelEtaReco == labelEtaMC ) // To be sure we compute the difference for the same particle
 499       {
 500         labelEtaReco = -1; // Restart of the label value to avoid double count the eta difference when computing the phi one
 501         Double_t EtaDif = EtaReco - EtaMC;
 502         MCHisto(eventSelection,triggerClassName,centrality,"EtaRes")->Fill(EtaDif);
 503         MCHisto(eventSelection,triggerClassName,centrality,"EtaResVsZ")->Fill(MCZv,EtaDif);
 504         MCHisto(eventSelection,triggerClassName,centrality,"EtaResVsnC")->Fill(nContributors,EtaDif);
 505       }
 506
 507          //Phi Resolution
 508       else if ( TString(p->GetName()).Contains("PhiReco") ) // We take the reco phi
 509       {
 510         sscanf(p->GetName(),"PhiReco%d",&labelPhiReco);
 511         PhiReco = p->GetVal();
 512         MCHisto(eventSelection,triggerClassName,centrality,"Phi")->Fill(PhiReco);
 513       }
 514       else if ( TString(p->GetName()).Contains("PhiMC") ) // We take the generated phi
 515       {
 516         sscanf(p->GetName(),"PhiMC%d",&labelPhiMC);
 517         PhiMC = p->GetVal();
 518         MCHisto(eventSelection,triggerClassName,centrality,"MCPhi")->Fill(PhiMC);
 519       }
 520
 521       if ( labelPhiReco > 0 && labelPhiReco == labelPhiMC ) // To be sure we compute the difference for the same particle
 522       {
 523         labelPhiReco = -1; // Restart of the label value to avoid double count the phi difference when computing the eta one
 524         Double_t PhiDif = PhiReco - PhiMC;
 525         MCHisto(eventSelection,triggerClassName,centrality,"PhiRes")->Fill(PhiDif);
 526
 527         //___With the following algorithm we refer the differences to the interval [-Pi/2,Pi/2]
 528         if ( PhiDif < -TMath::PiOver2() && PhiDif > -TMath::Pi() )
 529         {
 530           PhiDif = -TMath::Pi() - PhiDif;
 531         }
 532         else if ( PhiDif < -TMath::Pi() )
 533         {
 534           PhiDif = -2.*TMath::Pi() - PhiDif;
 535           if ( PhiDif < -TMath::PiOver2() )
 536           {
 537             PhiDif = -TMath::Pi() - PhiDif;
 538           }
 539         }
 540
 541         else if ( PhiDif > TMath::PiOver2() && PhiDif < TMath::Pi() )
 542         {
 543           PhiDif = TMath::Pi() - PhiDif;
 544         }
 545         else if ( PhiDif > TMath::Pi() )
 546         {
 547           PhiDif = 2.*TMath::Pi() - PhiDif;
 548           if ( PhiDif > TMath::PiOver2() )
 549           {
 550             PhiDif = TMath::Pi() - PhiDif;
 551           }
 552         }
 553         //___
 554
 555         MCHisto(eventSelection,triggerClassName,centrality,"PhiResVsZ")->Fill(MCZv,PhiDif);
 556         MCHisto(eventSelection,triggerClassName,centrality,"PhiResShifted")->Fill(PhiDif);
 557         MCHisto(eventSelection,triggerClassName,centrality,"PhiResVsnC")->Fill(nContributors,PhiDif);
 558       }
 559     }
 560
 561     return; // When computing resolutions we skip the rest of the method
 562   }
 563   //_____
 564
 565
 566   //___Multiplicity histos (just when not computing resolutions)___
 567   TH1* hSPDcorrectionVsEta = static_cast<TH1*>(nchList->FindObject("MCSPDcorrectionVsEta"));
 568   TH1* hNBkgTrackletsVSEta = static_cast<TH1*>(nchList->FindObject("NBkgTrackletsVSEta"));
 569   TH1* hNchVsPhi = static_cast<TH1*>(nchList->FindObject("MCNchVsPhi"));
 570   TH1* hNchVsEta = static_cast<TH1*>(nchList->FindObject("MCNchVsEta"));
 571   TH1* hNTrackletVsEta = static_cast<TH1*>(nchList->FindObject("MCNTrackletVsEta"));
 572   TH1* hNTrackletVsPhi = static_cast<TH1*>(nchList->FindObject("MCNTrackletVsPhi"));
 573
 574   TProfile* hMeanNchVsEta = static_cast<TProfile*>(MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
 575
 576   TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(MCHisto(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
 577
 578   Int_t nBins(0);
 579
 580   Double_t nchSum(0.),nTracklets(0.);
 581   for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) //Loop over eta bins
 582   {
 583     Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
 584
 585     if ( correction < 0 ) continue; // We count just the particles in the SPD acceptance.
 586
 587     nTracklets += hNTrackletVsEta->GetBinContent(j); // Reco tracklets
 588
 589     ++nBins; // We sum up the number of bins entering in the computation
 590
 591     Double_t eta = fEtaAxis->GetBinCenter(j);
 592     Double_t nch = hNchVsEta->GetBinContent(j); // Generated particles
 593
 594     nchSum += nch;
 595
 596     hMeanNchVsEta->Fill(eta,nch); // Fill the number of charged particles of each eta bin in the profile
 597
 598     hEventsVsZVertexVsEta->Fill(MCZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
 599   }
 600
 601   MCHisto(eventSelection,triggerClassName,centrality,"Tracklets")->Fill(nTracklets);
 602
 603   if ( vertex )
 604   {
 605     SPDZv = vertex->GetZ();
 606   }
 607
 608   Bool_t isMChisto = kTRUE; // Used to get the MC histos in the Add method
 609
 610   AddHisto(eventSelection,triggerClassName,centrality,"NBkgTrackletsVsZVertexVsEta",SPDZv,hNBkgTrackletsVSEta,isMChisto);
 611   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsPhi",MCZv,hNchVsPhi,isMChisto);
 612   AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",MCZv,hNchVsEta,isMChisto);
 613   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta",SPDZv,hNTrackletVsEta,isMChisto);
 614   AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi",SPDZv,hNTrackletVsPhi,isMChisto);
 615
 616   MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsZVertex")->Fill(MCZv,nchSum);
 617   MCHisto(eventSelection,triggerClassName,centrality,"Nch")->Fill(nchSum);
 618
 619
 620   // Mean dNch/dEta computation
 621   Double_t meandNchdEta(0.);
 622
 623   if ( nBins >  0 )
 624   {
 625     meandNchdEta = nchSum / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
 626   }
 627
 628   MCHisto(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
 629
 630   Double_t V0AMult = 0.;
 631   Double_t V0CMult = 0.;
 632
 633   AliVVZERO* vzero = Event()->GetVZEROData();
 634   if (vzero)
 635   {
 636     Double_t multV0A = vzero->GetMTotV0A();
 637     V0AMult = AliESDUtils::GetCorrV0A(multV0A,MCZv);
 638     Double_t multV0C = vzero->GetMTotV0C();
 639     V0CMult = AliESDUtils::GetCorrV0C(multV0C,MCZv);
 640
 641     MCHisto(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nchSum);
 642     MCHisto(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nchSum);
 643
 644   }
 645   //____
 646 }
 647
 648
 649 //_____________________________________________________________________________
 650 void AliAnalysisMuMuNch::GetEtaRangeSPD(Double_t spdZVertex, Double_t etaRange[])
 651 {
 652   // Returns the SPD eta range for a given z vertex.
 653
 654   Double_t etaMax(fEtaMax),etaMin(fEtaMin);
 655
 656   Double_t vf1LR = fSPD1LR->Eval(spdZVertex); //Eta values for the z vertex over the SPD acceptance curves
 657   Double_t vf2LR = fSPD2LR->Eval(spdZVertex);
 658   Double_t vf1LL = fSPD1LL->Eval(spdZVertex);
 659   Double_t vf2LL = fSPD2LL->Eval(spdZVertex);
 660
 661   //____We start by asigning the eta range as the eta values in the SPD acceptance curve
 662   if ( spdZVertex < 0. )
 663   {
 664     etaMax = vf2LR;
 665     etaMin = TMath::Max(vf1LL,vf2LL);
 666   }
 667   else
 668   {
 669     etaMax = TMath::Min(vf1LR,vf2LR);
 670     etaMin = vf2LL;
 671   }
 672   //____
 673
 674
 675   //____Algorithm to avoid taking bins which are crossed by the SPD acceptance curves
 676   Int_t binYMin = fEtaAxis->FindBin(etaMin); // Find the corresponding bins for eta max & min and z vertex
 677   Int_t binYMax = fEtaAxis->FindBin(etaMax);
 678   Int_t binX = fZAxis->FindBin(spdZVertex);
 679
 680   // Define the values for the relevant edges of the eta and z bins
 681   Double_t upEdge = fEtaAxis->GetBinUpEdge(binYMax);  // up edge of the top eta bin
 682   Double_t lowEdge = fEtaAxis->GetBinLowEdge(binYMin); // low edge of the bottom eta bin
 683   Double_t leftEdge = fZAxis->GetBinLowEdge(binX); // left edge of the z bin
 684   Double_t rightEdge = fZAxis->GetBinUpEdge(binX); // right edge of the z bin
 685
 686   Double_t etaMaxTemp(0.),etaMinTemp(0.);
 687   if ( spdZVertex < 0. )
 688   {
 689     etaMaxTemp = fSPD2LR->Eval(rightEdge); // Define the temporary eta max as the value of the curve int the righ edge of the bin
 690     etaMinTemp = TMath::Max(fSPD1LL->Eval(leftEdge),fSPD2LL->Eval(leftEdge));
 691   }
 692   else
 693   {
 694     etaMaxTemp = TMath::Min(fSPD1LR->Eval(rightEdge),fSPD2LR->Eval(rightEdge));
 695     etaMinTemp = fSPD2LL->Eval(leftEdge);
 696   }
 697
 698   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
 699   {
 700     binYMax = binYMax - 1; // Since the up edge of the current bin is bigger than the SPD acc curve we move 1 bin below
 701     upEdge = fEtaAxis->GetBinUpEdge(binYMax); // Take the up edge of the new bin
 702     etaMax = upEdge - 1E-6; // We substract 1E-6 cause the up edge of the a bin belongs to the bin+1
 703
 704   }
 705
 706   //We take eta min as the low edge of the 1st bin which is inside the SPD acceptance but not crossed by the curve
 707   while ( lowEdge < etaMinTemp )
 708   {
 709     binYMin = binYMin + 1; // Since the low edge of the current bin is smaller than the SPD acc curve we move 1 bin above
 710     lowEdge = fEtaAxis->GetBinLowEdge(binYMin); // Take the low edge of the new bin
 711     etaMin = lowEdge + 1E-6;
 712   }
 713   //____
 714
 715   // 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
 716   if ( etaMin < fEtaMin ) etaMin = fEtaMin + 1E-6;
 717   if ( etaMax > fEtaMax ) etaMax = fEtaMax - 1E-6;
 718
 719   etaRange[0] = etaMin;
 720   etaRange[1] = etaMax;
 721 }
 722
 723 //_____________________________________________________________________________
 724 Double_t AliAnalysisMuMuNch::GetSPDCorrection(Double_t zvert, Double_t eta) const
 725 {
 726   if (!fSPDCorrection)
 727   {
 728     AliFatal("ERROR: No SPD Correction");
 729     return 0;
 730   }
 731   Int_t bin = fSPDCorrection->FindBin(zvert,eta);
 732   return fSPDCorrection->GetBinContent(bin);
 733 }
 734
 735
 736 //==============================These methods are useless here, anyway they should go in the EventCutterClass=================
 737 //_____________________________________________________________________________
 738 Bool_t AliAnalysisMuMuNch::HasAtLeastNTrackletsInEtaRange(const AliVEvent& event, Int_t n, Double_t& etaMin, Double_t& etaMax) const
 739 {
 740   if ( event.IsA() != AliAODEvent::Class() )
 741   {
 742     return kFALSE;
 743   }
 744
 745   AliAODTracklets* tracklets = static_cast<const AliAODEvent&>(event).GetTracklets();
 746
 747   if (!tracklets)
 748   {
 749     return kFALSE;
 750   }
 751
 752   Int_t nTrackletsInRange(0);
 753
 754   Int_t nTracklets = tracklets->GetNumberOfTracklets();
 755
 756   for (Int_t i = 0 ; i < nTracklets && nTrackletsInRange < n; i++)
 757   {
 758     Double_t eta = -TMath::Log(TMath::Tan(tracklets->GetTheta(i)/2.0));
 759
 760     if ( eta > etaMin && eta < etaMax )
 761     {
 762       ++nTrackletsInRange;
 763     }
 764   }
 765
 766   return (nTrackletsInRange>=n);
 767 }
 768
 769 //_____________________________________________________________________________
 770 void AliAnalysisMuMuNch::NameOfHasAtLeastNTrackletsInEtaRange(TString& name, Int_t n, Double_t& etaMin, Double_t& etaMax) const
 771 {
 772   if ( TMath::AreEqualAbs(TMath::Abs(etaMin),TMath::Abs(etaMax),1E-9 ) )
 773   {
 774     name.Form("ATLEAST%dTRKLINABSETALT%3.1f",n,TMath::Abs(etaMin));
 775   }
 776   else
 777   {
 778     name.Form("ATLEAST%dTRKLINETA%3.1f-%3.1f",n,etaMin,etaMax);
 779   }
 780 }
 781
 782 //_____________________________________________________________________________
 783 void AliAnalysisMuMuNch::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
 784 {
 785   /// Set the event, compute multiplicities and add them to the event
 786
 787   if ( event->IsA() != AliAODEvent::Class() )
 788   {
 789     AliError("Only working for AODs for the moment.");
 790     return;
 791   }
 792
 793   AliAnalysisMuMuBase::SetEvent(event,mcEvent); // To have Event() and MCEvent() method working
 794
 795   TList* nchList = static_cast<TList*>(event->FindListObject("NCH")); // Define the list with the NCH info for each event
 796   if (!nchList)
 797   {
 798     nchList = new TList;
 799     nchList->SetOwner(kTRUE);
 800     nchList->SetName("NCH");
 801     event->AddObject(nchList);
 802   }
 803
 804   const AliAODVertex* vertexSPD = static_cast<const AliAODEvent*>(Event())->GetPrimaryVertexSPD(); // SPD vertex object
 805   AliAODTracklets* tracklets = static_cast<const AliAODEvent*>(Event())->GetTracklets(); // Tracklets object
 806
 807   TH1* hSPDcorrectionVsEta(0x0); // Pointers for the individual event histos
 808   TH1* hNchVsEta(0x0);
 809   TH1* hNTrackletVsEta(0x0);
 810   TH1* hNTrackletVsPhi(0x0);
 811
 812   //_______Create once the histos with the individual event "properties" (cleared at the beginning of each new event)_______
 813   if ( !fResolution ) // When computing resolutions we dont do anything else
 814   {
 815     if ( !Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta") )
 816     {
 817       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","SPDcorrectionVsEta","SPD correction-like vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 818
 819       CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsEta","Nch vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 820
 821       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsEta","Ntracklet vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 822
 823       CreateEventHistos(kHistoForData,"AliAnalysisMuMuNch","NBkgTrackletsVSEta","NBkg Tracklets vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 824
 825       Double_t phimin = 0.; //Phi range
 826       Double_t phimax = 2*TMath::Pi();
 827       Int_t nphibins = GetNbins(phimin,phimax,0.05);
 828
 829       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsPhi","Ntracklet vs #phi;#phi",nphibins,phimin,phimax);
 830
 831       CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsPhi","Nch vs #phi;#phi",nphibins,phimin,phimax);
 832
 833       CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","test","test",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
 834
 835       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD1LR);
 836       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD1LL);
 837       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD2LR);
 838       Histo("AliAnalysisMuMuNch","test")->GetListOfFunctions()->Add(fSPD2LL);
 839     }
 840     //_________
 841
 842
 843     hSPDcorrectionVsEta = Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta"); // Set the individual event histos
 844     hNTrackletVsEta = Histo("AliAnalysisMuMuNch","NTrackletVsEta");
 845     hNTrackletVsPhi = Histo("AliAnalysisMuMuNch","NTrackletVsPhi");
 846
 847
 848     hSPDcorrectionVsEta->Reset(); // Reset of the individual event histos
 849     hNTrackletVsEta->Reset();
 850     hNTrackletVsPhi->Reset();
 851   }
 852
 853   Double_t etaRange[2]; // Variables we will use in the multiplicity computation
 854   Int_t binMin,binMax;
 855   Int_t nTracklets(0);
 856   Double_t thetaTracklet(0.),phiTracklet(0.),etaTracklet(0.);
 857   Double_t nch(0.0);
 858   Int_t nBins(0);
 859
 860
 861   //____Data (or Reco if simu) multiplicity computation___
 862   if ( fSPDCorrection && !fResolution ) // When computing resolutions we dont do anything else
 863   {
 864     if ( tracklets && vertexSPD )
 865     {
 866       Double_t SPDZv = vertexSPD->GetZ();
 867
 868       GetEtaRangeSPD(SPDZv,etaRange);
 869
 870       Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,etaRange[1]);
 871       Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,etaRange[0]);
 872
 873       nTracklets = tracklets->GetNumberOfTracklets();
 874
 875       Double_t SPDr;
 876       for (Int_t i = 0 ; i < nTracklets ; i++)
 877       {
 878         thetaTracklet = tracklets->GetTheta(i);
 879         etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
 880         if ( etaTracklet < etaRange[0] || etaTracklet > etaRange[1] ) continue; // Avoid tracklets out of the eta SPD acceptance or out of the eta cut
 881
 882         SPDr = GetSPDCorrection(SPDZv,etaTracklet);
 883
 884         Int_t bin = fEtaAxis->FindBin(etaTracklet);
 885
 886         if ( SPDr!=0. && SPDr <= 2.5) // Threshold to reduce border effects
 887         {
 888           hSPDcorrectionVsEta->SetBinContent(bin,SPDr);
 889           hNTrackletVsEta->Fill(etaTracklet);
 890           hNTrackletVsPhi->Fill(tracklets->GetPhi(i));
 891         }
 892         else // If the correction is above the threshold we store a -1 in the correction to skip this eta bin int the fill method
 893         {
 894           hSPDcorrectionVsEta->SetBinContent(bin,-1.0);
 895         }
 896       }
 897
 898       //___ Fill the out-of-eta-range bins with -1.0
 899
 900       binMin = fEtaAxis->FindBin(etaRange[0]);
 901       binMax = fEtaAxis->FindBin(etaRange[1]);
 902
 903       for ( Int_t i = 1; i < binMin; ++i )
 904       {
 905         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
 906       }
 907       for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
 908       {
 909         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
 910       }
 911       //____
 912       nchList->Clear(); // We clear the NCH list for this new event
 913
 914       if ( !IsHistogrammingDisabled() )
 915       {
 916         nchList->Add(hSPDcorrectionVsEta->Clone());
 917         nchList->Add(hNTrackletVsEta->Clone());
 918         nchList->Add(hNTrackletVsPhi->Clone());
 919       }
 920
 921           //----- Mean dNchdEta computation to set it into the event
 922       for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) // Loop over eta bins
 923       {
 924         Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
 925
 926         Double_t eta = fEtaAxis->GetBinCenter(j);
 927
 928         if ( correction < 0 ) continue; // If the correction is < 0 we skip that eta bin
 929         else if ( correction == 0.0 ) // If the correction is 0 we have no tracklets in that eta bin
 930         {
 931           Double_t spdCorrection = GetSPDCorrection(SPDZv,eta);
 932           if ( spdCorrection == 0. || spdCorrection > 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)
 933         }
 934
 935         nch += hNTrackletVsEta->GetBinContent(j) * correction; // Number of charged particles (tracklets*SPDcorrection)
 936
 937         ++nBins; // We sum up the number of bins entering in the computation
 938       }
 939
 940       Double_t meandNchdEta(0.); // We compute the mean dNch/dEta in the event
 941       if ( nBins >  0 )
 942       {
 943         meandNchdEta = nch / (nBins*fEtaAxis->GetBinWidth(5));
 944       }
 945
 946       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.
 947       //------
 948
 949     }
 950
 951     else  nchList->Clear(); //To clear the NCH list for MC in case the event has no reconstructed SPD vertex
 952   }
 953   //_______
 954
 955   else nchList->Clear(); // To clear the NCH list for MC in case we have MC and no correction (SPD correction computation mode)
 956
 957
 958   //____Input MC multiplicity computation ___
 959   if ( HasMC() )
 960   {
 961     if ( !fResolution ) //When computing resolutions we dont do anything else
 962     {
 963       Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent());
 964       GetEtaRangeSPD(MCZv,etaRange);
 965
 966       hNchVsEta = MCHisto("AliAnalysisMuMuNch","NchVsEta");
 967       hSPDcorrectionVsEta = MCHisto("AliAnalysisMuMuNch","SPDcorrectionVsEta");
 968       TH1* hNchVsPhi = MCHisto("AliAnalysisMuMuNch","NchVsPhi");
 969       hNTrackletVsEta = MCHisto("AliAnalysisMuMuNch","NTrackletVsEta");
 970       hNTrackletVsPhi = MCHisto("AliAnalysisMuMuNch","NTrackletVsPhi");
 971
 972       hNTrackletVsEta->Reset();
 973       hNchVsEta->Reset();
 974       hNchVsPhi->Reset();
 975       hSPDcorrectionVsEta->Reset();
 976       hNTrackletVsPhi->Reset();
 977
 978       //___ Fill the out-of-eta-range bins with -1.0
 979       binMin = fEtaAxis->FindBin(etaRange[0]);
 980       binMax = fEtaAxis->FindBin(etaRange[1]);
 981
 982       for ( Int_t i = 1; i < binMin; ++i )
 983       {
 984         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
 985       }
 986       for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
 987       {
 988         hSPDcorrectionVsEta->SetBinContent(i,-1.0);
 989       }
 990       for ( Int_t i = binMin; i <= binMax; ++i ) // Fill the bins inside the eta range with +1
 991       {
 992         hSPDcorrectionVsEta->SetBinContent(i,1.0);
 993       }
 994       //___
 995
 996       Int_t nMCTracks = MCEvent()->GetNumberOfTracks(); // MC number of MC tracks
 997
 998       for ( Int_t i = 0; i < nMCTracks ; ++i ) //Loop over generated tracks
 999       {
1000         AliAODMCParticle* AODpart = static_cast<AliAODMCParticle*>(mcEvent->GetTrack(i));
1001
1002         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)
1003         {
1004           if ( AODpart->Charge()!=0 ) // We take only charged particles
1005           {
1006             hNchVsEta->Fill(AODpart->Eta());
1007             hNchVsPhi->Fill(AODpart->Phi());
1008
1009           }
1010         }
1011       }
1012
1013       nchList->Add(hNchVsEta->Clone("MCNchVsEta"));
1014       nchList->Add(hNchVsPhi->Clone("MCNchVsPhi"));
1015       nchList->Add(hSPDcorrectionVsEta->Clone("MCSPDcorrectionVsEta"));
1016     }
1017               //__Bkg tracklets and Resolution estimation __
1018     if ( tracklets ) // We can compute the Bkg tracklets and resolution only if we have the tracklets object in the event
1019     {
1020       TH1* hNBkgTrackletsVSEta(0x0); // Pointer for the Bkg histo
1021       if ( !fResolution )
1022       {
1023         hNBkgTrackletsVSEta = Histo("AliAnalysisMuMuNch","NBkgTrackletsVSEta");
1024
1025         hNBkgTrackletsVSEta->Reset();
1026       }
1027
1028       nTracklets = tracklets->GetNumberOfTracklets();
1029
1030       for (Int_t i = 0 ; i < nTracklets ; i++) // Loop over tracklets to check if they come or not from the same MC particle
1031       {
1032         thetaTracklet = tracklets->GetTheta(i);
1033         etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
1034         phiTracklet = tracklets->GetPhi(i);
1035
1036
1037         if ( !fResolution )
1038         {
1039           hNTrackletVsEta->Fill(etaTracklet);
1040           hNTrackletVsPhi->Fill(phiTracklet);
1041         }
1042
1043         Int_t label1 = tracklets->GetLabel(i,0);
1044         Int_t label2 = tracklets->GetLabel(i,1);
1045
1046         if (label1 != label2 ) // Tracklets not comming from the same MC particle are Bkg
1047         {
1048          if (!fResolution ) hNBkgTrackletsVSEta->Fill(etaTracklet);
1049         }
1050         else if ( !fSPDCorrection && fResolution ) // Compute the resolutions with the tracklets comming from the same MC particle
1051         {
1052           AliAODMCParticle* AODpartMC = static_cast<AliAODMCParticle*>(MCEvent()->GetTrack(label1));
1053           Double_t etaTrackletMC = AODpartMC->Eta();
1054           Double_t phiTrackletMC = AODpartMC->Phi();
1055
1056           // Resolution variables
1057           nchList->Add(new TParameter<Double_t>(Form("EtaReco%d",label1),etaTracklet));
1058           nchList->Add(new TParameter<Double_t>(Form("EtaMC%d",label1),etaTrackletMC));
1059
1060           nchList->Add(new TParameter<Double_t>(Form("PhiReco%d",label1),phiTracklet));
1061           nchList->Add(new TParameter<Double_t>(Form("PhiMC%d",label1),phiTrackletMC));
1062
1063         }
1064       }
1065
1066       if (!fResolution )
1067       {
1068         nchList->Add(hNTrackletVsEta->Clone("MCNTrackletVsEta"));
1069         nchList->Add(hNTrackletVsPhi->Clone("MCNTrackletVsPhi"));
1070         nchList->Add(hNBkgTrackletsVSEta->Clone());
1071       }
1072
1073     }
1074
1075   }
1076   //_______
1077
1078 }
1079
1080 //_____________________________________________________________________________
1081 void AliAnalysisMuMuNch::Terminate(Option_t *)
1082 {
1083   /// Called once at the end of the query
1084   if ( !HistogramCollection() ) return;
1085
1086   if ( HistogramCollection()->FindObject("/MCINPUT/AliAnalysisMuMuNch/NTrackletVsEta") )
1087   {
1088     HistogramCollection()->Remove("/MCINPUT/AliAnalysisMuMuNch/NTrackletVsEta");
1089     HistogramCollection()->Remove("/MCINPUT/AliAnalysisMuMuNch/NTrackletVsPhi");
1090     HistogramCollection()->Remove("/MCINPUT/AliAnalysisMuMuNch/NchVsEta");
1091     HistogramCollection()->Remove("/MCINPUT/AliAnalysisMuMuNch/NchVsPhi");
1092     HistogramCollection()->Remove("/MCINPUT/AliAnalysisMuMuNch/SPDcorrectionVsEta");
1093     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NBkgTrackletsVSEta");
1094   }
1095
1096   if ( HistogramCollection()->FindObject("/AliAnalysisMuMuNch/NTrackletVsEta") )
1097   {
1098     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsEta");
1099     HistogramCollection()->Remove("/AliAnalysisMuMuNch/test");
1100     HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsPhi");
1101     HistogramCollection()->Remove("/AliAnalysisMuMuNch/SPDcorrectionVsEta");
1102   }
1103   //____ Compute dNchdEta histo
1104   TObjArray* idArr =  HistogramCollection()->SortAllIdentifiers();
1105
1106   TIter next(idArr);
1107   TObjString* id;
1108
1109   while ( (id = static_cast<TObjString*>(next())) )
1110   {
1111     TProfile* p = static_cast<TProfile*>(HistogramCollection()->FindObject(Form("%s%s",id->GetName(),"MeanNchVsEta")));
1112
1113     if ( !p ) continue;
1114
1115     TH1* h = new TH1F("MeandNchdEta","Event averaged dN_{ch}/d#eta ;#eta;<dN_{ch}/d#eta>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1116
1117     if ( p->GetNbinsX() != h->GetNbinsX() || p->GetXaxis()->GetXmin() != h->GetXaxis()->GetXmin() || p->GetXaxis()->GetXmax() != h->GetXaxis()->GetXmax() )
1118     {
1119       AliError("ERROR: Cannot compute MeandNchdEta since the binning doesn't match with MeanNchVsEta histo");
1120       continue;
1121     }
1122
1123     for ( Int_t i = 1 ; i < h->GetNbinsX() ; i++ )
1124     {
1125       h->SetBinContent(i,p->GetBinContent(i)/p->GetBinWidth(i));
1126       h->SetBinError(i,p->GetBinError(i)/p->GetBinWidth(i));
1127       h->SetEntries(p->GetEntries());
1128     }
1129
1130     HistogramCollection()->Adopt(Form("%s",id->GetName()),h);
1131   }
1132
1133   delete idArr;
1134   //___
1135 }