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