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