1 #include "AliAnalysisMuMuNch.h"
5 * \ingroup pwg-muon-mumu
7 * \class AliAnalysisMuMuNch
9 * SPD tracklet to Nch analysis
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.
17 #include "AliAODTracklets.h"
18 #include "AliAnalysisMuonUtility.h"
20 #include "AliAODEvent.h"
21 #include "AliESDUtils.h"
23 #include "AliAnalysisMuMuCutRegistry.h"
24 #include "AliAnalysisMuMuCutElement.h"
25 #include "Riostream.h"
26 #include "TParameter.h"
29 #include "AliMergeableCollection.h"
30 #include "AliMCEvent.h"
31 #include "AliAODMCParticle.h"
37 #include "TObjString.h"
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
44 // par[0] = radius of SPD layer
53 theta = TMath::ATan(par[0]/z);
56 else if( d > 14.09999 )
59 theta = TMath::Pi() - TMath::ATan(par[0]/z);
62 return -TMath::Log(TMath::Tan(theta/2.));
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
67 // par[0] = radius of SPD layer
76 theta = TMath::Pi() - TMath::ATan(par[0]/z);
82 theta = TMath::ATan(par[0]/z);
85 return -TMath::Log(TMath::Tan(theta/2.));
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(),
95 fEtaAxis(new TAxis(TMath::Nint(10./0.1),-5.,5.)),
96 fZAxis(new TAxis(TMath::Nint(80/0.25),-40.,40.)),
102 fResolution(computeResolution),
110 fSPDCorrection = static_cast<TH2F*>(spdCorrection->Clone());
111 fSPDCorrection->SetDirectory(0);
113 DefineSPDAcceptance();
117 DisableHistograms("*");
121 //_____________________________________________________________________________
122 AliAnalysisMuMuNch::~AliAnalysisMuMuNch()
124 delete fSPDCorrection;
133 //_____________________________________________________________________________
134 void AliAnalysisMuMuNch::DefineHistogramCollection(const char* eventSelection,
135 const char* triggerClassName,
136 const char* centrality)
138 // Define multiplicity histos
140 if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuNch") )
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);
148 Double_t multMin = 0.; //Tracklets multiplicity range
149 Double_t multMax = 500.;
150 Int_t nbinsMult = GetNbins(multMin,multMax,1.);
152 Double_t phimin = 0.; //Phi range
153 Double_t phimax = 2*TMath::Pi();
154 Int_t nphibins = GetNbins(phimin,phimax,0.05);
156 if ( !fSPDCorrection && fResolution ) // Resolution histos
158 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"EtaRes","#eta resolution;#eta_{Reco} - #eta_{MC};Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
160 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiRes","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*2,-phimax,phimax);
162 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"PhiResShifted","#phi resolution;#phi_{Reco} - #phi_{MC};Counts",nphibins*4,-phimax/4,phimax/4);
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());
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);
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());
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);
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());
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());
176 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Phi","Reco #phi distribution;#phi;Counts",nphibins,phimin,phimax);
178 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCPhi","MC #phi distribution;#phi;Counts",nphibins,phimin,phimax);
180 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","Reco #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
182 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"MCEta","MC #eta distribution;#phi;Counts",(fEtaAxis->GetNbins())*2,fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
184 return; // When computing resolutions we don't want to create the rest of histos
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
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");
193 CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Tracklets","Number of tracklets distribution;N_{Tracklets};N_{events}",nbinsMult,multMin,multMax);
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");
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");
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);
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");
206 CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Nch","Number of charged particles distribution;N_{ch};N_{events}",nbinsMult,multMin,multMax);
208 CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"dNchdEta","<dNchdEta> distribution;dN_{ch}/d#eta;N_{events}",nbinsMult,multMin,multMax);
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
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);
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);
216 // profile histograms
218 CreateEventHistos(kHistoForData,eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta","Mean number of tracklets vs #eta;#eta;<N_{Tracklets}>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax(),0);
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);
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)
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);
228 //_____________________________________________________________________________
229 void AliAnalysisMuMuNch::DefineSPDAcceptance()
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);
244 //_____________________________________________________________________________
245 void AliAnalysisMuMuNch::AddHisto(const char* eventSelection,
246 const char* triggerClassName,
247 const char* centrality,
248 const char* histoname,
252 // Adds the content of a 1D histo to a 2D histo a the z position
254 Int_t zbin = fZAxis->FindBin(z);
256 if (isMC) h2 = static_cast<TH2F*>(MCHisto(eventSelection,triggerClassName,centrality,histoname));
257 else h2 = static_cast<TH2F*>(Histo(eventSelection,triggerClassName,centrality,histoname));
259 for ( Int_t i = 1; i <= h->GetXaxis()->GetNbins(); ++i )
261 Double_t content = h2->GetCellContent(zbin,i);
263 if ( h->GetBinContent(i) > 0 )
265 h2->SetCellContent(zbin,i,content + h->GetBinContent(i));
269 h2->SetEntries(h2->GetSumOfWeights());
273 //_____________________________________________________________________________
274 void AliAnalysisMuMuNch::AttachSPDAcceptance(UInt_t dataType,
275 const char* eventSelection,
276 const char* triggerClassName,
277 const char* centrality,const char* histoname)
279 if ( dataType & kHistoForData )
281 if( !Histo(eventSelection,triggerClassName,centrality,histoname) )
283 AliError(Form("ERROR: SPD Acceptance curves attach failed. Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
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);
292 if ( dataType & kHistoForMCInput )
294 if( !MCHisto(eventSelection,triggerClassName,centrality,histoname) )
296 AliError(Form("ERROR: SPD Acceptance curves attach failed. MC Histo /%s/%s/%s/%s not found",eventSelection,triggerClassName,centrality,histoname));
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);
309 //_____________________________________________________________________________
310 void AliAnalysisMuMuNch::FillHistosForEvent(const char* eventSelection,
311 const char* triggerClassName,
312 const char* centrality)
314 // Fills the data (or reco if simu) multiplicity histos
316 if ( IsHistogrammingDisabled() ) return;
318 if ( fResolution ) return; //When computing resolutions we skip this method
320 if ( !AliAnalysisMuonUtility::IsAODEvent(Event()) )
322 AliError("Don't know how to deal with ESDs...");
326 if ( HasMC() && !fSPDCorrection ) // We have MC but no correction (SPD correction computation mode)so we skip the method
331 AliAODEvent* aod = static_cast<AliAODEvent*>(Event());
333 AliVVertex* vertex = aod->GetPrimaryVertexSPD();
335 TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
337 if (!nchList || nchList->IsEmpty() ) // Empty NCH means that there is no SPD vertex ( see SetEvent() ) when runing on data.
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
349 else SPDZv = vertex->GetZ();
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"));
355 TH1* hNchVsEta =static_cast<TH1*>(hNTrackletVsEta->Clone("NchVsEta"));
357 TProfile* hMeanTrackletsVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanTrackletsVsEta"));
358 TProfile* hMeanNchVsEta = static_cast<TProfile*>(Histo(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
360 TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(Histo(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
365 Double_t nTracklets(0.0);
367 for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) // Loop over eta bins
369 Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
371 Double_t eta = fEtaAxis->GetBinCenter(j);
373 if ( correction < 0 ) continue;
374 else if ( correction == 0.0 ) // No tracklets found in this eta bin.
376 correction = GetSPDCorrection(SPDZv,eta);
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
381 Double_t ntr = hNTrackletVsEta->GetBinContent(j); // Tracklets in eta bin
383 nch += ntr * correction; // Number of charged particles (corrected tracklets)
387 hMeanTrackletsVsEta->Fill(eta,ntr); // Fill the number of tracklets of each eta bin in the profile
389 hMeanNchVsEta->Fill(eta,ntr*correction);
391 ++nBins; // We sum up the number of bins entering in the computation
393 // Fill the number of charged particles of each eta bin in the profile
394 hNchVsEta->SetBinContent( j, hNchVsEta->GetBinContent(j) * correction );
396 hEventsVsZVertexVsEta->Fill(SPDZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
399 AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsPhi",SPDZv,hNTrackletVsPhi);
400 AddHisto(eventSelection,triggerClassName,centrality,"TrackletsVsZVertexVsEta",SPDZv,hNTrackletVsEta);
402 AddHisto(eventSelection,triggerClassName,centrality,"NchVsZVertexVsEta",SPDZv,hNchVsEta);
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);
408 Histo(eventSelection,triggerClassName,centrality,"TrackletsVsNch")->Fill(nTracklets,nch);
409 Histo(eventSelection,triggerClassName,centrality,"Nch")->Fill(nch);
411 Double_t V0AMult = 0.;
412 Double_t V0CMult = 0.;
414 AliVVZERO* vzero = aod->GetVZEROData();
417 Double_t multV0A = vzero->GetMTotV0A();
418 V0AMult = AliESDUtils::GetCorrV0A(multV0A,SPDZv);
419 Double_t multV0C = vzero->GetMTotV0C();
420 V0CMult = AliESDUtils::GetCorrV0C(multV0C,SPDZv);
422 Histo(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nch);
423 Histo(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nch);
430 // Mean dNch/dEta computation
431 Double_t meandNchdEta(0.);
435 meandNchdEta = nch / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
438 Histo(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
444 //_____________________________________________________________________________
445 void AliAnalysisMuMuNch::FillHistosForMCEvent(const char* eventSelection,const char* triggerClassName,const char* centrality)
447 /// Fill input MC multiplicity histos
449 if ( IsHistogrammingDisabled() ) return;
451 TList* nchList = static_cast<TList*>(Event()->FindListObject("NCH"));
453 if (!nchList || nchList->IsEmpty())
458 Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent()); // Definition of MC generated z vertex
459 AliVVertex* vertex = static_cast<AliAODEvent*>(Event())->GetPrimaryVertexSPD();
463 //____Resolution Histos___
464 if ( !fSPDCorrection && fResolution )
466 Int_t nContributors(0);
469 nContributors = vertex->GetNContributors();
470 SPDZv = vertex->GetZ();
472 MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsnC")->Fill(nContributors,SPDZv - MCZv);
473 MCHisto(eventSelection,triggerClassName,centrality,"SPDZvResVsMCz")->Fill(MCZv,SPDZv - MCZv);
475 Double_t EtaReco(0.),EtaMC(0.),PhiReco(0.),PhiMC(0.);
476 Int_t i(-1),labelEtaReco(-1),labelEtaMC(-1),labelPhiReco(-1),labelPhiMC(-1);
478 while ( i < nchList->GetEntries() - 1 )
481 while ( nchList->At(i)->IsA() != TObjString::Class() ) // In case there is a diferent object, just to skip it
486 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(nchList->At(i));
490 if ( TString(p->GetName()).Contains("EtaReco") ) // We take the reco eta
492 sscanf(p->GetName(),"EtaReco%d",&labelEtaReco);
493 EtaReco = p->GetVal();
494 MCHisto(eventSelection,triggerClassName,centrality,"Eta")->Fill(EtaReco);
496 else if ( TString(p->GetName()).Contains("EtaMC") ) // We take the generated eta
498 sscanf(p->GetName(),"EtaMC%d",&labelEtaMC);
500 MCHisto(eventSelection,triggerClassName,centrality,"MCEta")->Fill(EtaMC);
502 if ( labelEtaReco > 0 && labelEtaReco == labelEtaMC ) // To be sure we compute the difference for the same particle
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);
512 else if ( TString(p->GetName()).Contains("PhiReco") ) // We take the reco phi
514 sscanf(p->GetName(),"PhiReco%d",&labelPhiReco);
515 PhiReco = p->GetVal();
516 MCHisto(eventSelection,triggerClassName,centrality,"Phi")->Fill(PhiReco);
518 else if ( TString(p->GetName()).Contains("PhiMC") ) // We take the generated phi
520 sscanf(p->GetName(),"PhiMC%d",&labelPhiMC);
522 MCHisto(eventSelection,triggerClassName,centrality,"MCPhi")->Fill(PhiMC);
525 if ( labelPhiReco > 0 && labelPhiReco == labelPhiMC ) // To be sure we compute the difference for the same particle
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);
531 //___With the following algorithm we refer the differences to the interval [-Pi/2,Pi/2]
532 if ( PhiDif < -TMath::PiOver2() && PhiDif > -TMath::Pi() )
534 PhiDif = -TMath::Pi() - PhiDif;
536 else if ( PhiDif < -TMath::Pi() )
538 PhiDif = -2.*TMath::Pi() - PhiDif;
539 if ( PhiDif < -TMath::PiOver2() )
541 PhiDif = -TMath::Pi() - PhiDif;
545 else if ( PhiDif > TMath::PiOver2() && PhiDif < TMath::Pi() )
547 PhiDif = TMath::Pi() - PhiDif;
549 else if ( PhiDif > TMath::Pi() )
551 PhiDif = 2.*TMath::Pi() - PhiDif;
552 if ( PhiDif > TMath::PiOver2() )
554 PhiDif = TMath::Pi() - PhiDif;
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);
565 return; // When computing resolutions we skip the rest of the method
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"));
578 TProfile* hMeanNchVsEta = static_cast<TProfile*>(MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsEta"));
580 TH2* hEventsVsZVertexVsEta = static_cast<TH2*>(MCHisto(eventSelection,triggerClassName,centrality,"EventsVsZVertexVsEta"));
584 Double_t nchSum(0.),nTracklets(0.);
585 for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) //Loop over eta bins
587 Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
589 if ( correction < 0 ) continue; // We count just the particles in the SPD acceptance.
591 nTracklets += hNTrackletVsEta->GetBinContent(j); // Reco tracklets
593 ++nBins; // We sum up the number of bins entering in the computation
595 Double_t eta = fEtaAxis->GetBinCenter(j);
596 Double_t nch = hNchVsEta->GetBinContent(j); // Generated particles
600 hMeanNchVsEta->Fill(eta,nch); // Fill the number of charged particles of each eta bin in the profile
602 hEventsVsZVertexVsEta->Fill(MCZv,eta,1.0); // Fill 1 count each eta bin where the events contributes
605 MCHisto(eventSelection,triggerClassName,centrality,"Tracklets")->Fill(nTracklets);
609 SPDZv = vertex->GetZ();
612 Bool_t isMChisto = kTRUE; // Used to get the MC histos in the Add method
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);
620 MCHisto(eventSelection,triggerClassName,centrality,"MeanNchVsZVertex")->Fill(MCZv,nchSum);
621 MCHisto(eventSelection,triggerClassName,centrality,"Nch")->Fill(nchSum);
624 // Mean dNch/dEta computation
625 Double_t meandNchdEta(0.);
629 meandNchdEta = nchSum / (nBins*fEtaAxis->GetBinWidth(5)); // Divide by nBins to get the mean and by the binWidht to get the d/dEta
632 MCHisto(eventSelection,triggerClassName,centrality,"dNchdEta")->Fill(meandNchdEta);
634 Double_t V0AMult = 0.;
635 Double_t V0CMult = 0.;
637 AliVVZERO* vzero = Event()->GetVZEROData();
640 Double_t multV0A = vzero->GetMTotV0A();
641 V0AMult = AliESDUtils::GetCorrV0A(multV0A,MCZv);
642 Double_t multV0C = vzero->GetMTotV0C();
643 V0CMult = AliESDUtils::GetCorrV0C(multV0C,MCZv);
645 MCHisto(eventSelection,triggerClassName,centrality,"V0AMultVsNch")->Fill(V0AMult,nchSum);
646 MCHisto(eventSelection,triggerClassName,centrality,"V0CMultVsNch")->Fill(V0CMult,nchSum);
653 //_____________________________________________________________________________
654 void AliAnalysisMuMuNch::GetEtaRangeSPD(Double_t spdZVertex, Double_t etaRange[])
656 // Returns the SPD eta range for a given z vertex.
658 Double_t etaMax(fEtaMax),etaMin(fEtaMin);
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);
665 //____We start by asigning the eta range as the eta values in the SPD acceptance curve
666 if ( spdZVertex < 0. )
669 etaMin = TMath::Max(vf1LL,vf2LL);
673 etaMax = TMath::Min(vf1LR,vf2LR);
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);
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
690 Double_t etaMaxTemp(0.),etaMinTemp(0.);
691 if ( spdZVertex < 0. )
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));
698 etaMaxTemp = TMath::Min(fSPD1LR->Eval(rightEdge),fSPD2LR->Eval(rightEdge));
699 etaMinTemp = fSPD2LL->Eval(leftEdge);
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
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
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 )
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;
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;
723 etaRange[0] = etaMin;
724 etaRange[1] = etaMax;
727 //_____________________________________________________________________________
728 Double_t AliAnalysisMuMuNch::GetSPDCorrection(Double_t zvert, Double_t eta) const
732 AliFatal("ERROR: No SPD Correction");
735 Int_t bin = fSPDCorrection->FindBin(zvert,eta);
736 return fSPDCorrection->GetBinContent(bin);
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
744 if ( event.IsA() != AliAODEvent::Class() )
749 AliAODTracklets* tracklets = static_cast<const AliAODEvent&>(event).GetTracklets();
756 Int_t nTrackletsInRange(0);
758 Int_t nTracklets = tracklets->GetNumberOfTracklets();
760 for (Int_t i = 0 ; i < nTracklets && nTrackletsInRange < n; i++)
762 Double_t eta = -TMath::Log(TMath::Tan(tracklets->GetTheta(i)/2.0));
764 if ( eta > etaMin && eta < etaMax )
770 return (nTrackletsInRange>=n);
773 //_____________________________________________________________________________
774 void AliAnalysisMuMuNch::NameOfHasAtLeastNTrackletsInEtaRange(TString& name, Int_t n, Double_t& etaMin, Double_t& etaMax) const
776 if ( TMath::AreEqualAbs(TMath::Abs(etaMin),TMath::Abs(etaMax),1E-9 ) )
778 name.Form("ATLEAST%dTRKLINABSETALT%3.1f",n,TMath::Abs(etaMin));
782 name.Form("ATLEAST%dTRKLINETA%3.1f-%3.1f",n,etaMin,etaMax);
786 //_____________________________________________________________________________
787 void AliAnalysisMuMuNch::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
789 /// Set the event, compute multiplicities and add them to the event
791 if ( event->IsA() != AliAODEvent::Class() )
793 AliError("Only working for AODs for the moment.");
797 AliAnalysisMuMuBase::SetEvent(event,mcEvent); // To have Event() and MCEvent() method working
799 TList* nchList = static_cast<TList*>(event->FindListObject("NCH")); // Define the list with the NCH info for each event
803 nchList->SetOwner(kTRUE);
804 nchList->SetName("NCH");
805 event->AddObject(nchList);
808 const AliAODVertex* vertexSPD = static_cast<const AliAODEvent*>(Event())->GetPrimaryVertexSPD(); // SPD vertex object
809 AliAODTracklets* tracklets = static_cast<const AliAODEvent*>(Event())->GetTracklets(); // Tracklets object
811 TH1* hSPDcorrectionVsEta(0x0); // Pointers for the individual event histos
813 TH1* hNTrackletVsEta(0x0);
814 TH1* hNTrackletVsPhi(0x0);
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
819 if ( !Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta") )
821 CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","SPDcorrectionVsEta","SPD correction-like vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
823 CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsEta","Nch vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
825 CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsEta","Ntracklet vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
827 CreateEventHistos(kHistoForData,"AliAnalysisMuMuNch","NBkgTrackletsVSEta","NBkg Tracklets vs #eta;#eta",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
829 Double_t phimin = 0.; //Phi range
830 Double_t phimax = 2*TMath::Pi();
831 Int_t nphibins = GetNbins(phimin,phimax,0.05);
833 CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","NTrackletVsPhi","Ntracklet vs #phi;#phi",nphibins,phimin,phimax);
835 CreateEventHistos(kHistoForMCInput,"AliAnalysisMuMuNch","NchVsPhi","Nch vs #phi;#phi",nphibins,phimin,phimax);
837 CreateEventHistos(kHistoForData | kHistoForMCInput,"AliAnalysisMuMuNch","test","test",fZAxis->GetNbins(),fZAxis->GetXmin(),fZAxis->GetXmax(),fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
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);
847 hSPDcorrectionVsEta = Histo("AliAnalysisMuMuNch","SPDcorrectionVsEta"); // Set the individual event histos
848 hNTrackletVsEta = Histo("AliAnalysisMuMuNch","NTrackletVsEta");
849 hNTrackletVsPhi = Histo("AliAnalysisMuMuNch","NTrackletVsPhi");
852 hSPDcorrectionVsEta->Reset(); // Reset of the individual event histos
853 hNTrackletVsEta->Reset();
854 hNTrackletVsPhi->Reset();
857 Double_t etaRange[2]; // Variables we will use in the multiplicity computation
860 Double_t thetaTracklet(0.),phiTracklet(0.),etaTracklet(0.);
865 //____Data (or Reco if simu) multiplicity computation___
866 if ( fSPDCorrection && !fResolution ) // When computing resolutions we dont do anything else
868 if ( tracklets && vertexSPD )
870 Double_t SPDZv = vertexSPD->GetZ();
872 GetEtaRangeSPD(SPDZv,etaRange);
874 Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,etaRange[1]);
875 Histo("AliAnalysisMuMuNch","test")->Fill(SPDZv,etaRange[0]);
877 nTracklets = tracklets->GetNumberOfTracklets();
880 for (Int_t i = 0 ; i < nTracklets ; i++)
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
886 SPDr = GetSPDCorrection(SPDZv,etaTracklet);
888 Int_t bin = fEtaAxis->FindBin(etaTracklet);
890 if ( SPDr!=0. && SPDr <= 2.5) // Threshold to reduce border effects
892 hSPDcorrectionVsEta->SetBinContent(bin,SPDr);
893 hNTrackletVsEta->Fill(etaTracklet);
894 hNTrackletVsPhi->Fill(tracklets->GetPhi(i));
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
898 hSPDcorrectionVsEta->SetBinContent(bin,-1.0);
902 //___ Fill the out-of-eta-range bins with -1.0
904 binMin = fEtaAxis->FindBin(etaRange[0]);
905 binMax = fEtaAxis->FindBin(etaRange[1]);
907 for ( Int_t i = 1; i < binMin; ++i )
909 hSPDcorrectionVsEta->SetBinContent(i,-1.0);
911 for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
913 hSPDcorrectionVsEta->SetBinContent(i,-1.0);
916 nchList->Clear(); // We clear the NCH list for this new event
918 if ( !IsHistogrammingDisabled() )
920 nchList->Add(hSPDcorrectionVsEta->Clone());
921 nchList->Add(hNTrackletVsEta->Clone());
922 nchList->Add(hNTrackletVsPhi->Clone());
925 //----- Mean dNchdEta computation to set it into the event
926 for (Int_t j = 1 ; j <= fEtaAxis->GetNbins() ; j++) // Loop over eta bins
928 Double_t correction = hSPDcorrectionVsEta->GetBinContent(j);
930 Double_t eta = fEtaAxis->GetBinCenter(j);
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
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)
939 nch += hNTrackletVsEta->GetBinContent(j) * correction; // Number of charged particles (tracklets*SPDcorrection)
941 ++nBins; // We sum up the number of bins entering in the computation
944 Double_t meandNchdEta(0.); // We compute the mean dNch/dEta in the event
947 meandNchdEta = nch / (nBins*fEtaAxis->GetBinWidth(5));
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.
955 else nchList->Clear(); //To clear the NCH list for MC in case the event has no reconstructed SPD vertex
959 else nchList->Clear(); // To clear the NCH list for MC in case we have MC and no correction (SPD correction computation mode)
962 //____Input MC multiplicity computation ___
965 if ( !fResolution ) //When computing resolutions we dont do anything else
967 Double_t MCZv = AliAnalysisMuonUtility::GetMCVertexZ(Event(),MCEvent());
968 GetEtaRangeSPD(MCZv,etaRange);
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");
976 hNTrackletVsEta->Reset();
979 hSPDcorrectionVsEta->Reset();
980 hNTrackletVsPhi->Reset();
982 //___ Fill the out-of-eta-range bins with -1.0
983 binMin = fEtaAxis->FindBin(etaRange[0]);
984 binMax = fEtaAxis->FindBin(etaRange[1]);
986 for ( Int_t i = 1; i < binMin; ++i )
988 hSPDcorrectionVsEta->SetBinContent(i,-1.0);
990 for ( Int_t i = binMax + 1 ; i <= fEtaAxis->GetNbins(); ++i )
992 hSPDcorrectionVsEta->SetBinContent(i,-1.0);
994 for ( Int_t i = binMin; i <= binMax; ++i ) // Fill the bins inside the eta range with +1
996 hSPDcorrectionVsEta->SetBinContent(i,1.0);
1000 Int_t nMCTracks = MCEvent()->GetNumberOfTracks(); // MC number of MC tracks
1002 for ( Int_t i = 0; i < nMCTracks ; ++i ) //Loop over generated tracks
1004 AliAODMCParticle* AODpart = static_cast<AliAODMCParticle*>(mcEvent->GetTrack(i));
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)
1008 if ( AODpart->Charge()!=0 ) // We take only charged particles
1010 hNchVsEta->Fill(AODpart->Eta());
1011 hNchVsPhi->Fill(AODpart->Phi());
1017 nchList->Add(hNchVsEta->Clone("MCNchVsEta"));
1018 nchList->Add(hNchVsPhi->Clone("MCNchVsPhi"));
1019 nchList->Add(hSPDcorrectionVsEta->Clone("MCSPDcorrectionVsEta"));
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
1024 TH1* hNBkgTrackletsVSEta(0x0); // Pointer for the Bkg histo
1027 hNBkgTrackletsVSEta = Histo("AliAnalysisMuMuNch","NBkgTrackletsVSEta");
1029 hNBkgTrackletsVSEta->Reset();
1032 nTracklets = tracklets->GetNumberOfTracklets();
1034 for (Int_t i = 0 ; i < nTracklets ; i++) // Loop over tracklets to check if they come or not from the same MC particle
1036 thetaTracklet = tracklets->GetTheta(i);
1037 etaTracklet = -TMath::Log(TMath::Tan(thetaTracklet/2.));
1038 phiTracklet = tracklets->GetPhi(i);
1043 hNTrackletVsEta->Fill(etaTracklet);
1044 hNTrackletVsPhi->Fill(phiTracklet);
1047 Int_t label1 = tracklets->GetLabel(i,0);
1048 Int_t label2 = tracklets->GetLabel(i,1);
1050 if (label1 != label2 ) // Tracklets not comming from the same MC particle are Bkg
1052 if (!fResolution ) hNBkgTrackletsVSEta->Fill(etaTracklet);
1054 else if ( !fSPDCorrection && fResolution ) // Compute the resolutions with the tracklets comming from the same MC particle
1056 AliAODMCParticle* AODpartMC = static_cast<AliAODMCParticle*>(MCEvent()->GetTrack(label1));
1057 Double_t etaTrackletMC = AODpartMC->Eta();
1058 Double_t phiTrackletMC = AODpartMC->Phi();
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));
1064 nchList->Add(new TParameter<Double_t>(Form("PhiReco%d",label1),phiTracklet));
1065 nchList->Add(new TParameter<Double_t>(Form("PhiMC%d",label1),phiTrackletMC));
1072 nchList->Add(hNTrackletVsEta->Clone("MCNTrackletVsEta"));
1073 nchList->Add(hNTrackletVsPhi->Clone("MCNTrackletVsPhi"));
1074 nchList->Add(hNBkgTrackletsVSEta->Clone());
1084 //_____________________________________________________________________________
1085 void AliAnalysisMuMuNch::Terminate(Option_t *)
1087 /// Called once at the end of the query
1088 if ( !HistogramCollection() ) return;
1090 if ( HistogramCollection()->FindObject("/MCINPUT/AliAnalysisMuMuNch/NTrackletVsEta") )
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");
1100 if ( HistogramCollection()->FindObject("/AliAnalysisMuMuNch/NTrackletVsEta") )
1102 HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsEta");
1103 HistogramCollection()->Remove("/AliAnalysisMuMuNch/test");
1104 HistogramCollection()->Remove("/AliAnalysisMuMuNch/NTrackletVsPhi");
1105 HistogramCollection()->Remove("/AliAnalysisMuMuNch/SPDcorrectionVsEta");
1107 //____ Compute dNchdEta histo
1108 TObjArray* idArr = HistogramCollection()->SortAllIdentifiers();
1113 while ( (id = static_cast<TObjString*>(next())) )
1115 TProfile* p = static_cast<TProfile*>(HistogramCollection()->FindObject(Form("%s%s",id->GetName(),"MeanNchVsEta")));
1119 TH1* h = new TH1F("MeandNchdEta","Event averaged dN_{ch}/d#eta ;#eta;<dN_{ch}/d#eta>",fEtaAxis->GetNbins(),fEtaAxis->GetXmin(),fEtaAxis->GetXmax());
1121 if ( p->GetNbinsX() != h->GetNbinsX() || p->GetXaxis()->GetXmin() != h->GetXaxis()->GetXmin() || p->GetXaxis()->GetXmax() != h->GetXaxis()->GetXmax() )
1123 AliError("ERROR: Cannot compute MeandNchdEta since the binning doesn't match with MeanNchVsEta histo");
1127 for ( Int_t i = 1 ; i < h->GetNbinsX() ; i++ )
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());
1134 HistogramCollection()->Adopt(Form("%s",id->GetName()),h);