]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muon/AliAnalysisMuMuMinv.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliAnalysisMuMu classes update
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisMuMuMinv.cxx
CommitLineData
5376e016
CP		 1#include "AliAnalysisMuMuMinv.h"
2
3/**
4 *
5 * \ingroup pwg-muon-mumu
6 *
7 * \class AliAnalysisMuMuMinv
16560e8e 8 *
9 * Analysis which fills a bunch of histograms for invariant mass analysis of J/psi
5376e016
CP		 10 *
11 * Can be used on real data and/or simulated (for instance to get Acc x Eff)
12 *
13 * Can optionally use as input an already computed Acc x Eff matrix that will be applied
14 * when filling the invariant mass histograms.
15 *
16 */
17
18#include "TH2F.h"
19#include "AliLog.h"
20#include "TObjArray.h"
21#include "AliAnalysisMuMuBinning.h"
22#include "TString.h"
23#include "TLorentzVector.h"
24#include "TString.h"
25#include "AliMCEvent.h"
26#include "AliMergeableCollection.h"
27#include "AliAnalysisMuonUtility.h"
28#include "TParameter.h"
29#include <cassert>
30
31ClassImp(AliAnalysisMuMuMinv)
32
33//_____________________________________________________________________________
16560e8e 34AliAnalysisMuMuMinv::AliAnalysisMuMuMinv(TH2* accEffHisto, Bool_t computeMeanPt, Int_t systLevel)
5376e016 35: AliAnalysisMuMuBase(),
16560e8e 36fcomputeMeanPt(computeMeanPt),
37fAccEffHisto(0x0),
38fsystLevel(systLevel)
5376e016
CP		 39{
40 // FIXME : find the AccxEff histogram from HistogramCollection()->Histo("/EXCHANGE/JpsiAccEff")
41
16560e8e 42 if ( accEffHisto )
43 {
44 fAccEffHisto = static_cast<TH2F*>(accEffHisto->Clone());
45 fAccEffHisto->SetDirectory(0);
46 }
5376e016
CP		 47}
48
49//_____________________________________________________________________________
50AliAnalysisMuMuMinv::~AliAnalysisMuMuMinv()
51{
52 /// dtor
53 delete fAccEffHisto;
54}
55
56//_____________________________________________________________________________
57void
58AliAnalysisMuMuMinv::DefineHistogramCollection(const char* eventSelection,
59 const char* triggerClassName,
60 const char* centrality)
61{
62 /// Define the histograms this analysis will use
63
64 if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv") )
65 {
66 return;
67 }
68
69 // dummy histogram to signal that we already defined all our histograms (see above)
70 CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv","Dummy semaphore",1,0,1);
71
72 /// Create invariant mass histograms
73
74 Double_t minvMin = 0;
75 Double_t minvMax = 16;
76 Int_t nMinvBins = GetNbins(minvMin,minvMax,0.025);
77
78 Int_t nMCMinvBins = GetNbins(minvMin,minvMax,0.1);
79
16560e8e 80 Double_t rapidityMin = -5;
81 Double_t rapidityMax = -2;
82 Int_t nbinsRapidity = GetNbins(rapidityMin,rapidityMax,0.05);
5376e016 83
16560e8e 84 Double_t etaMin = -5;
85 Double_t etaMax = -2;
86 Int_t nbinsEta = GetNbins(etaMin,etaMax,0.05);
5376e016 87
16560e8e 88 TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr","");//We may include ,v0a,v0acent
89
90 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Pt","#mu+#mu- Pt distribution",
91 200,0,20,-2);
92
93 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Y","#mu+#mu- Y distribution",
94 nbinsRapidity,rapidityMin,rapidityMax,-2);
95
96 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","#mu+#mu- Eta distribution",
97 nbinsEta,etaMin,etaMax);
98
99
100 //___Histos for pure MC
101 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Pt","MCINPUT #mu+#mu- Pt distribution",
102 200,0,20,-2);
103
104 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Y","MCINPUT #mu+#mu- Y distribution",
105 nbinsRapidity,rapidityMin,rapidityMax,-2);
106
107 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","MCINPUT #mu+#mu- Eta distribution",
108 nbinsEta,etaMin,etaMax);
109
110 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Pt","MCINPUT #mu+#mu- Pt distribution",
111 200,0,20,-2);
112
113 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Y","MCINPUT #mu+#mu- Y distribution",
114 nbinsRapidity,rapidityMin,rapidityMax,-2);
115
116 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Eta","MCINPUT #mu+#mu- Eta distribution",
117 nbinsEta,etaMin,etaMax);
118 //____
119
5376e016
CP		 120 // CreatePairHistos(eventSelection,triggerClassName,centrality,"BinFlowPt","#mu+#mu- BinFlowPt distribution",
121 // 200,0,20);
122
16560e8e 123 CreatePairHistos(kHistoForData,eventSelection,triggerClassName,centrality,"PtRecVsSim","#mu+#mu- Pt distribution rec vs sim",
5376e016
CP		 124 200,0,20,200,0,20);
125
5376e016
CP		 126 TIter next(bins);
127 AliAnalysisMuMuBinning::Range* r;
128 Int_t nb(0);
129
130 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
131 {
16560e8e 132 TString minvName(GetMinvHistoName(*r,kFALSE));
5376e016
CP		 133
134 ++nb;
135
16560e8e 136 if ( !IsHistogramDisabled(minvName.Data()) )
5376e016 137 {
5376e016 138
16560e8e 139 AliDebug(1,Form("bin %d %s histoname = %s",nb,r->AsString().Data(),minvName.Data()));
140
141 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
142 Form("#mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
143 r->AsString().Data()),
144 nMinvBins,minvMin,minvMax,-2);
145
146 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
147 Form("MCINPUT #mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
148 r->AsString().Data()),
149 nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
150
151 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),minvName.Data(),
152 Form("MCINPUT #mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
153 r->AsString().Data()),
154 nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
155
156
157 if ( fcomputeMeanPt )
158 {
159 TString mPtName(Form("MeanPtVs%s",minvName.Data()));
160 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
161 Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
162 r->AsString().Data()),
163 nMinvBins,minvMin,minvMax,0);
164
165 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
166 Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
167 r->AsString().Data()),
168 nMinvBins,minvMin,minvMax,0); //Pure MC Histo
169
170 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),mPtName.Data(),
171 Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
172 r->AsString().Data()),
173 nMinvBins,minvMin,minvMax,0); //Pure MC Histo
174 }
175
176// if ( HasMC() )
177// {
178// TH1* h = new TH1F(minvName.Data(),Form("MC #mu+#mu- inv. mass %s",r->AsString().Data()),
179// nMCMinvBins,minvMin,minvMax);
180//
181// HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),h);
182//
183// HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),static_cast<TH1*>(h->Clone()));
184// }
185 }
186
187 if ( ShouldCorrectDimuonForAccEff() )
188 {
189 minvName = GetMinvHistoName(*r,kTRUE);
5376e016 190
16560e8e 191 if ( !IsHistogramDisabled(minvName.Data()) )
192 {
193
194 AliDebug(1,Form("bin %d %s histoname = %s",nb,r->AsString().Data(),minvName.Data()));
195
196 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
197 Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}}(GeV/c^2);Counts",
198 r->AsString().Data()),
199 nMinvBins,minvMin,minvMax,-2);
200
201 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
202 Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
203 r->AsString().Data()),
204 nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
205
206 CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),minvName.Data(),
207 Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
208 r->AsString().Data()),
209 nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
210
211 if ( fcomputeMeanPt )
212 {
213 TString mPtName(Form("MeanPtVs%s",minvName.Data()));
214 CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
215 Form("#mu+#mu- mean p_{T} %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-}}>",
216 r->AsString().Data()),
217 nMinvBins,minvMin,minvMax,0);
218 }
219
220// if ( HasMC() )
221// {
222// TH1* h = new TH1F(minvName.Data(),Form("MC #mu+#mu- inv. mass %s",r->AsString().Data()),
223// nMCMinvBins,minvMin,minvMax);
224//
225// h->Sumw2();
226//
227// HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),h);
228//
229// HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),static_cast<TH1*>(h->Clone()));
230// }
231 }
5376e016
CP		 232 }
233 }
16560e8e 234
5376e016
CP		 235 delete bins;
236}
237
238//_____________________________________________________________________________
239void AliAnalysisMuMuMinv::FillHistosForPair(const char* eventSelection, const char* triggerClassName,
240 const char* centrality, const char* pairCutName,
241 const AliVParticle& tracki,
242 const AliVParticle& trackj)
243{
5376e016 244
5376e016
CP		 245 TLorentzVector pi(tracki.Px(),tracki.Py(),tracki.Pz(),
246 TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+tracki.P()*tracki.P()));
247
16560e8e 248 if (!AliAnalysisMuonUtility::IsMuonTrack(&tracki) ) return;
249 if (!AliAnalysisMuonUtility::IsMuonTrack(&trackj) ) return;
5376e016
CP		 250
251 TLorentzVector pair4Momentum(trackj.Px(),trackj.Py(),trackj.Pz(),
252 TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+trackj.P()*trackj.P()));
253
254 pair4Momentum += pi;
255
256
16560e8e 257//// if (!IsHistogramDisabled("Chi12"))
258//// {
259//// Histo(eventSelection,triggerClassName,centrality,pairCutName,"Chi12")
260//// ->Fill(
261//// AliAnalysisMuonUtility::GetChi2perNDFtracker(&tracki),
262//// AliAnalysisMuonUtility::GetChi2perNDFtracker(&trackj));
263//// }
264////
265//// if (!IsHistogramDisabled("Rabs12"))
266//// {
267//// Histo(eventSelection,triggerClassName,centrality,pairCutName,"Rabs12")
268//// ->Fill(AliAnalysisMuonUtility::GetRabs(&tracki),
269//// AliAnalysisMuonUtility::GetRabs(&trackj));
270//// }
5376e016
CP		 271
272 if ( ( tracki.Charge() != trackj.Charge() ) )
273 {
16560e8e 274 Double_t inputWeight = WeightDistribution(pair4Momentum.Pt(),pair4Momentum.Rapidity());
275
5376e016
CP		 276 if ( !IsHistogramDisabled("Pt") )
277 {
16560e8e 278 Histo(eventSelection,triggerClassName,centrality,pairCutName,"Pt")->Fill(pair4Momentum.Pt(),inputWeight);
5376e016 279 }
16560e8e 280 if ( !IsHistogramDisabled("Y") )
5376e016 281 {
16560e8e 282 Histo(eventSelection,triggerClassName,centrality,pairCutName,"Y")->Fill(pair4Momentum.Rapidity(),inputWeight);
5376e016 283 }
16560e8e 284 if ( !IsHistogramDisabled("Eta") )
285 {
286 Histo(eventSelection,triggerClassName,centrality,pairCutName,"Eta")->Fill(pair4Momentum.Eta());
287 }
288
289 TLorentzVector* pair4MomentumMC(0x0);
290 Double_t inputWeightMC(1.);
5376e016
CP		 291 if ( HasMC() )
292 {
293 Int_t labeli = tracki.GetLabel();
294 Int_t labelj = trackj.GetLabel();
295
296 if ( labeli < 0 || labelj < 0 )
297 {
298 AliError("Got negative labels!");
299 }
300 else
301 {
302 AliVParticle* mcTracki = MCEvent()->GetTrack(labeli);
303 AliVParticle* mcTrackj = MCEvent()->GetTrack(labelj);
304
305 TLorentzVector mcpi(mcTracki->Px(),mcTracki->Py(),mcTracki->Pz(),
306 TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTracki->P()*mcTracki->P()));
307 TLorentzVector mcpj(mcTrackj->Px(),mcTrackj->Py(),mcTrackj->Pz(),
308 TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTrackj->P()*mcTrackj->P()));
309
310 mcpj += mcpi;
311
16560e8e 312 inputWeightMC = WeightDistribution(mcpj.Pt(),mcpj.Rapidity());
313
5376e016
CP		 314 Histo(eventSelection,triggerClassName,centrality,pairCutName,"PtRecVsSim")->Fill(mcpj.Pt(),pair4Momentum.Pt());
315
16560e8e 316 if ( !IsHistogramDisabled("Pt") )
317 {
318 MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Pt")->Fill(mcpj.Pt(),inputWeightMC);
319 }
320 if ( !IsHistogramDisabled("Y") )
321 {
322 MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Y")->Fill(mcpj.Rapidity(),inputWeightMC);
323 }
324 if ( !IsHistogramDisabled("Eta") )
325 {
326 MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Eta")->Fill(mcpj.Eta());
327 }
328 pair4MomentumMC = &mcpj;
329
5376e016
CP		 330 }
331 }
332
16560e8e 333
334 TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr",""); // We may include: ,v0a,v0acent
5376e016
CP		 335 TIter nextBin(bins);
336 AliAnalysisMuMuBinning::Range* r;
337
338 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
339 {
340 Bool_t ok(kFALSE);
16560e8e 341 Bool_t okMC(kFALSE);
5376e016
CP		 342
343 if ( r->IsIntegrated() )
344 {
345 ok = kTRUE;
16560e8e 346 if ( pair4MomentumMC ) okMC = kTRUE;
5376e016
CP		 347 }
348 else if ( r->Is2D() )
349 {
350 if ( r->AsString().BeginsWith("PTVSY") )
351 {
352 ok = r->IsInRange(pair4Momentum.Rapidity(),pair4Momentum.Pt());
16560e8e 353 if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Rapidity(),pair4MomentumMC->Pt());
5376e016
CP		 354 }
355 else if ( r->AsString().BeginsWith("YVSPT") )
356 {
357 ok = r->IsInRange(pair4Momentum.Pt(),pair4Momentum.Rapidity());
16560e8e 358 if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity());
359 }
360 else if ( r->Quantity() == "NTRCORRPT" )
361 {
362 TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
363 if (list)
364 {
365 Int_t i(-1);
366 Bool_t parFound(kFALSE);
367 while ( i < list->GetEntries() - 1 && !parFound )
368 {
369 i++;
370 while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
371 {
372 i++;
373 }
374
375 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
376
377 if ( TString(p->GetName()).Contains("NtrCorr") )
378 {
379 parFound = kTRUE;
380 ok = r->IsInRange(p->GetVal(),pair4Momentum.Pt());
381 }
382 }
383 }
384
385 }
386 else if ( r->Quantity() == "NTRCORRY" )
387 {
388 TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
389 if (list)
390 {
391 Int_t i(-1);
392 Bool_t parFound(kFALSE);
393 while ( i < list->GetEntries() - 1 && !parFound )
394 {
395 i++;
396 while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
397 {
398 i++;
399 }
400
401 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
402
403 if ( TString(p->GetName()).Contains("NtrCorr") )
404 {
405 parFound = kTRUE;
406 ok = r->IsInRange(p->GetVal(),pair4Momentum.Rapidity());
407 }
408 }
409 }
410
5376e016 411 }
16560e8e 412
5376e016
CP		 413 else
414 {
415 AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
416 }
417 }
418 else
419 {
420 if ( r->Quantity() == "PT" )
421 {
422 ok = r->IsInRange(pair4Momentum.Pt());
16560e8e 423 if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Pt());
5376e016
CP		 424 }
425 else if ( r->Quantity() == "Y" )
426 {
427 ok = r->IsInRange(pair4Momentum.Rapidity());
16560e8e 428 if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Rapidity());
5376e016
CP		 429 }
430 else if ( r->Quantity() == "PHI" )
431 {
432 ok = r->IsInRange(pair4Momentum.Phi());
16560e8e 433 if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Phi());
5376e016 434 }
16560e8e 435 else if ( r->Quantity() == "DNCHDETA" )
5376e016
CP		 436 {
437 TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
438 if (list)
439 {
16560e8e 440 Int_t i(-1);
441 Bool_t parFound(kFALSE);
442 while ( i < list->GetEntries() - 1 && !parFound )
5376e016 443 {
16560e8e 444 i++;
445 while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
446 {
447 i++;
448 }
449
450 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
451
452 if ( TString(p->GetName()).Contains("dNchdEta") )
5376e016 453 {
16560e8e 454 parFound = kTRUE;
5376e016 455 ok = r->IsInRange(p->GetVal());
5376e016
CP		 456 }
457 }
458 }
16560e8e 459// else AliFatal("No dNchdEta info on Event");
460
461 }
462 else if ( r->Quantity() == "NTRCORR" || r->Quantity() == "RELNTRCORR" )
463 {
464 TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
465 if (list)
466 {
467 Int_t i(-1);
468 Bool_t parFound(kFALSE);
469 while ( i < list->GetEntries() - 1 && !parFound )
470 {
471 i++;
472 while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
473 {
474 i++;
475 }
476
477 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
478
479 if ( TString(p->GetName()).Contains("NtrCorr") )
480 {
481 parFound = kTRUE;
482 if ( r->Quantity() == "NTRCORR" ) ok = r->IsInRange(p->GetVal());
483 else ok = r->IsInRange(p->GetVal()/5.97);
484 }
485 }
486 }
487// else AliFatal("No ntrcorr info on Event");
488
5376e016
CP		 489 }
490 }
491
16560e8e 492 if ( ok || okMC )
5376e016
CP		 493 {
494 TString minvName = GetMinvHistoName(*r,kFALSE);
495
496 if (!IsHistogramDisabled(minvName.Data()))
497 {
16560e8e 498 TH1* h(0x0);
499 if ( ok )
500 {
501 h = Histo(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
502
503 if (!h)
504 {
505 AliError(Form("Could not get %s",minvName.Data()));
506 //continue;
507 }
508 else h->Fill(pair4Momentum.M(),inputWeight);
509 }
510 if( okMC )
511 {
512 h = MCHisto(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
513
514 if (!h)
515 {
516 AliError(Form("Could not get MC %s",minvName.Data()));
517 //continue;
518 }
519 else h->Fill(pair4MomentumMC->M(),inputWeightMC);
520 }
521
522 if ( fcomputeMeanPt )
523 {
524 TString hprofName("");
525
526 if ( ok )
527 {
528 hprofName= Form("MeanPtVs%s",minvName.Data());
529
530 TProfile* hprof = Prof(eventSelection,triggerClassName,centrality,pairCutName,hprofName.Data());
531
532 if ( !hprof )
533 {
534 AliError(Form("Could not get %s",hprofName.Data()));
535 }
536
537 else
538 {
539 // hprof->Approximate(); //I dont think its necessary here
540 hprof->Fill(pair4Momentum.M(),pair4Momentum.Pt(),inputWeight);
541 }
542 }
543 if ( okMC )
544 {
545 hprofName= Form("MeanPtVs%s",minvName.Data());
546
547 TProfile* hprof = MCProf(eventSelection,triggerClassName,centrality,pairCutName,hprofName.Data());
548
549 if ( !hprof )
550 {
551 AliError(Form("Could not get MC %s",hprofName.Data()));
552 }
553
554 else
555 {
556 // hprof->Approximate(); //I dont think its necessary here
557 hprof->Fill(pair4MomentumMC->M(),pair4MomentumMC->Pt(),inputWeightMC);
558 }
559 }
560
561 }
562
563 }
564
565 if ( ShouldCorrectDimuonForAccEff() )
566 {
567
568 Double_t AccxEff(0);
569 Bool_t okAccEff(kFALSE);
570 if ( ok )
571 {
572 AccxEff = GetAccxEff(pair4Momentum.Pt(),pair4Momentum.Rapidity());
573 if ( AccxEff <= 0.0 )
574 {
575 AliError(Form("AccxEff < 0 for pt = %f & y = %f ",pair4Momentum.Pt(),pair4Momentum.Rapidity()));
576 // continue;
577 }
578 else okAccEff = kTRUE;
579 }
5376e016 580
16560e8e 581 Double_t AccxEffMC(0);
582 Bool_t okAccEffMC(kFALSE);
583 if ( okMC )
5376e016 584 {
16560e8e 585 AccxEffMC= GetAccxEff(pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity());
586 if ( AccxEffMC <= 0.0 )
587 {
588 AliError(Form("AccxEff < 0 for MC pair with pt = %f & y = %f ",pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity()));
589 // continue;
590 }
591 else okAccEffMC = kTRUE;
5376e016 592 }
5376e016 593
16560e8e 594 minvName = GetMinvHistoName(*r,kTRUE);
595
596 if (!IsHistogramDisabled(minvName.Data()))
5376e016 597 {
16560e8e 598 TH1* hCorr = Histo(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
599
600 if (!hCorr)
601 {
602 AliError(Form("Could not get %sr",minvName.Data()));
603 }
604
605 else if ( okAccEff ) hCorr->Fill(pair4Momentum.M(),inputWeight/AccxEff);
606
607 if( okAccEffMC )
608 {
609 hCorr = MCHisto(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
610
611 if (!hCorr)
612 {
613 AliError(Form("Could not get MC %s",minvName.Data()));
614 //continue;
615 }
616 else hCorr->Fill(pair4MomentumMC->M(),inputWeightMC/AccxEffMC);
617 }
618
619 if ( fcomputeMeanPt )
620 {
621 TString hprofCorrName("");
622 if( ok )
623 {
624 hprofCorrName = Form("MeanPtVs%s",minvName.Data());
625
626 TProfile* hprofCorr = Prof(eventSelection,triggerClassName,centrality,pairCutName,hprofCorrName.Data());
627
628 if ( !hprofCorr )
629 {
630 AliError(Form("Could not get %s",hprofCorrName.Data()));
631 }
632 else if ( okAccEff )
633 {
634 // hprofCorr->Approximate(); //I dont know if its necessary here
635 hprofCorr->Fill(pair4Momentum.M(),pair4Momentum.Pt(),inputWeight/AccxEff);
636 }
637 }
638 if( okMC )
639 {
640 hprofCorrName = Form("MeanPtVs%s",minvName.Data());
641
642 TProfile* hprofCorr = MCProf(eventSelection,triggerClassName,centrality,pairCutName,hprofCorrName.Data());
643
644 if ( !hprofCorr )
645 {
646 AliError(Form("Could not get MC %s",hprofCorrName.Data()));
647 }
648 else if ( okAccEffMC )
649 {
650 // hprofCorr->Approximate(); //I dont know if its necessary here
651 hprofCorr->Fill(pair4MomentumMC->M(),pair4MomentumMC->Pt(),inputWeightMC/AccxEffMC);
652 }
653 }
654
655 }
656
5376e016 657 }
16560e8e 658
5376e016
CP		 659 }
660 }
661 }
662
663 delete bins;
664 }
16560e8e 665
5376e016
CP		 666}
667
16560e8e 668
5376e016 669//_____________________________________________________________________________
16560e8e 670void AliAnalysisMuMuMinv::FillHistosForMCEvent(const char* eventSelection,const char* triggerClassName,const char* centrality)
5376e016 671{
16560e8e 672 // Fill the input Monte-Carlo histograms related to muons. Intended to be used on pure simulations so we wont use eventSelection, triggerClassName and centrality variables.
673
674 if ( !HasMC() ) return;
675
676 // Specific things for MC // These histos should go in the AliAnalysisMuMuGlobal task, but then we have to loop 2 times on the MCTracks...
677// if (!Histo(MCInputPrefix(),"ALL","Pt"))
678// {
679// HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Pt","Pt",200,0,25));
680// HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Pt","Pt",200,0,25));
681//
682// Double_t rapidityMin = -5;
683// Double_t rapidityMax = -2;
684// Int_t nbinsRapidity = GetNbins(rapidityMin,rapidityMax,0.05);
685//
686// HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
687// HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
688//
689// Double_t etaMin = -5;
690// Double_t etaMax = -2;
691// Int_t nbinsEta = GetNbins(etaMin,etaMax,0.05);
692//
693// HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
694// HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
695// }
5376e016
CP		 696
697 Int_t nMCTracks = MCEvent()->GetNumberOfTracks();
16560e8e 698
699 TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr","");//We may include: ,v0a,v0acent.
700 TIter nextBin(bins);
701 AliAnalysisMuMuBinning::Range* r;
702
5376e016
CP		 703 for ( Int_t i = 0; i < nMCTracks; ++i )
704 {
705 AliVParticle* part = MCEvent()->GetTrack(i);
706
16560e8e 707 if (AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
708 AliAnalysisMuonUtility::GetMotherIndex(part)==-1)
5376e016 709 {
16560e8e 710 Double_t inputWeight = WeightDistribution(part->Pt(),part->Y());
711
712 MCHisto(eventSelection,triggerClassName,centrality,"Pt")->Fill(part->Pt(),inputWeight);
713 MCHisto(eventSelection,triggerClassName,centrality,"Y")->Fill(part->Y(),inputWeight);
714 MCHisto(eventSelection,triggerClassName,centrality,"Eta")->Fill(part->Eta());
715
716 if ( part->Y() < -2.5 && part->Y() > -4.0 )
717 {
718 MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Pt")->Fill(part->Pt(),inputWeight);
719 MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Y")->Fill(part->Y(),inputWeight);
720 MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Eta")->Fill(part->Eta());
721 }
722
723 nextBin.Reset();
724
725 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
726 {
727 Bool_t ok(kFALSE);
728
729 if ( r->IsIntegrated() )
730 {
731 ok = kTRUE;
732 }
733 else if ( r->Is2D() )
734 {
735 if ( r->AsString().BeginsWith("PTVSY") )
736 {
737 ok = r->IsInRange(part->Y(),part->Pt());
738 }
739 else if ( r->AsString().BeginsWith("YVSPT") )
740 {
741 ok = r->IsInRange(part->Pt(),part->Y());
742 }
743 else
744 {
745 AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
746 }
747 }
748 else
749 {
750 if ( r->Quantity() == "PT" )
751 {
752 ok = r->IsInRange(part->Pt());
753 }
754 else if ( r->Quantity() == "Y" )
755 {
756 ok = r->IsInRange(part->Y());
757 }
758 else if ( r->Quantity() == "PHI" )
759 {
760 ok = r->IsInRange(part->Phi());
761 }
762 }
763
764 if ( ok )
765 {
766 TString hname = GetMinvHistoName(*r,kFALSE);
767
768 if (!IsHistogramDisabled(hname.Data()))
769 {
770
771 TH1* h = MCHisto(eventSelection,triggerClassName,centrality,hname.Data());
772
773 if (!h)
774 {
775 AliError(Form("Could not get /%s/%s/%s/%s/ %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
776 continue;
777 }
778
779 h->Fill(part->M(),inputWeight);
780
781 if ( part->Y() < -2.5 && part->Y() > -4.0 )
782 {
783 h = MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data());
784 if (!h)
785 {
786 AliError(Form("Could not get /%s/%s/%s/%s/INYRANGE %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
787 continue;
788 }
789 h->Fill(part->M(),inputWeight);
790 }
791
792 }
793
794 if ( fcomputeMeanPt )
795 {
796 TString hprofName= Form("MeanPtVs%s",hname.Data());
797
798 TProfile* hprof = MCProf(eventSelection,triggerClassName,centrality,hprofName.Data());
799
800 if ( !hprof )
801 {
802 AliError(Form("Could not get %s",hprofName.Data()));
803 }
804 else
805 {
806 // hprof->Approximate(); //I dont think its necessary here
807 hprof->Fill(part->M(),part->Pt(),inputWeight);
808 }
809
5376e016
CP		 810 if ( part->Y() < -2.5 && part->Y() > -4.0 )
811 {
16560e8e 812 hprof = MCProf(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hprofName.Data());
813
814 if ( !hprof )
815 {
816 AliError(Form("Could not get %s",hprofName.Data()));
817 }
818 else
819 {
820 // hprof->Approximate(); //I dont think its necessary here
821 hprof->Fill(part->M(),part->Pt(),inputWeight);
822 }
823
5376e016
CP		 824 }
825
16560e8e 826
827 }
828
829// if ( ShouldCorrectDimuonForAccEff() ) //What is the sense of this?? We should not correct the input
830// {
831// Double_t AccxEff = GetAccxEff(part->Pt(),part->Y());
832// if ( AccxEff <= 0.0 )
833// {
834// AliError(Form("AccxEff < 0 for pt = %f & y = %f ",part->Pt(),part->Y()));
835// continue;
836// }
837// hname = GetMinvHistoName(*r,kTRUE);
838//
839// if (!IsHistogramDisabled(hname.Data()))
840// {
841//
842// TH1* h = MCHisto(eventSelection,triggerClassName,centrality,hname.Data());
843//
844// if (!h)
845// {
846// AliError(Form("Could not get /%s/%s/%s/%s/ %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
847// continue;
848// }
849//
850// h->Fill(part->M(),inputWeight/AccxEff);
851//
852// if ( part->Y() < -2.5 && part->Y() > -4.0 )
853// {
854// h = MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data());
855// if (!h)
856// {
857// AliError(Form("Could not get /%s/%s/%s/%s/INYRANGE %s",MCInputPrefix(),eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data()));
858// continue;
859// }
860// h->Fill(part->M(),inputWeight./AccxEff);
861// }
862//
863// }
864//
865// }
866
867 }
868 }
5376e016 869 }
5376e016
CP		 870 }
871
16560e8e 872 delete bins;
5376e016 873
16560e8e 874
5376e016
CP		 875// for ( Int_t i = 0; i < nMCTracks; ++i )
876// {
877// AliVParticle* part = MCEvent()->GetTrack(i);
878//
5376e016
CP		 879// if (AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
880// AliAnalysisMuonUtility::GetMotherIndex(part)==-1)
881// {
882//
16560e8e 883// Histo(MCInputPrefix(),"ALL","Pt")->Fill(part->Pt());
884// Histo(MCInputPrefix(),"ALL","Y")->Fill(part->Y());
885// Histo(MCInputPrefix(),"ALL","Eta")->Fill(part->Eta());
5376e016
CP		 886//
887// if ( part->Y() < -2.5 && part->Y() > -4.0 )
888// {
16560e8e 889// Histo(MCInputPrefix(),"INYRANGE","Pt")->Fill(part->Pt());
890// Histo(MCInputPrefix(),"INYRANGE","Y")->Fill(part->Y());
891// Histo(MCInputPrefix(),"INYRANGE","Eta")->Fill(part->Eta());
5376e016
CP		 892// }
893//
894// nextBin.Reset();
895//
896// while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
897// {
898// Bool_t ok(kFALSE);
899//
900// if ( r->IsIntegrated() )
901// {
902// ok = kTRUE;
903// }
904// else if ( r->Is2D() )
905// {
906// if ( r->AsString().BeginsWith("PTVSY") )
907// {
908// ok = r->IsInRange(part->Y(),part->Pt());
909// }
910// else if ( r->AsString().BeginsWith("YVSPT") )
911// {
912// ok = r->IsInRange(part->Pt(),part->Y());
913// }
914// else
915// {
916// AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
917// }
918// }
919// else
920// {
921// if ( r->Quantity() == "PT" )
922// {
923// ok = r->IsInRange(part->Pt());
924// }
925// else if ( r->Quantity() == "Y" )
926// {
927// ok = r->IsInRange(part->Y());
928// }
929// else if ( r->Quantity() == "PHI" )
930// {
931// ok = r->IsInRange(part->Phi());
932// }
933// }
934//
935// if ( ok )
936// {
937// TString hname = GetMinvHistoName(*r,kFALSE);
938//
939// if (!IsHistogramDisabled(hname.Data()))
940// {
941//
16560e8e 942// TH1* h = Histo(MCInputPrefix(),"ALL",hname.Data());
5376e016
CP		 943//
944// if (!h)
945// {
946// AliError(Form("Could not get ALL %s",hname.Data()));
947// continue;
948// }
949//
950// h->Fill(part->M());
951//
952// if ( part->Y() < -2.5 && part->Y() > -4.0 )
953// {
16560e8e 954// h = Histo(MCInputPrefix(),"INYRANGE",hname.Data());
5376e016
CP		 955// if (!h)
956// {
957// AliError(Form("Could not get INYRANGE %s",hname.Data()));
958// continue;
959// }
960// h->Fill(part->M());
961// }
962//
963// }
964//
16560e8e 965// if ( ShouldCorrectDimuonForAccEff() )
966// {
967// Double_t AccxEff = GetAccxEff(part->Pt(),part->Y());
968// if ( AccxEff <= 0.0 )
969// {
970// AliError(Form("AccxEff < 0 for pt = %f & y = %f ",part->Pt(),part->Y()));
971// continue;
972// }
973// hname = GetMinvHistoName(*r,kTRUE);
974//
975// if (!IsHistogramDisabled(hname.Data()))
976// {
977//
978// TH1* h = Histo(MCInputPrefix(),"ALL",hname.Data());
979//
980// if (!h)
981// {
982// AliError(Form("Could not get ALL %s",hname.Data()));
983// continue;
984// }
985//
986// h->Fill(part->M(),1./AccxEff);
987//
988// if ( part->Y() < -2.5 && part->Y() > -4.0 )
989// {
990// h = Histo(MCInputPrefix(),"INYRANGE",hname.Data());
991// if (!h)
992// {
993// AliError(Form("Could not get INYRANGE %s",hname.Data()));
994// continue;
995// }
996// h->Fill(part->M(),1./AccxEff);
997// }
998//
999// }
1000//
1001// }
1002//
5376e016
CP		 1003// }
1004// }
1005// }
16560e8e 1006//
1007// delete bins;
5376e016 1008// }
16560e8e 1009}
5376e016
CP		 1010
1011//_____________________________________________________________________________
1012TString AliAnalysisMuMuMinv::GetMinvHistoName(const AliAnalysisMuMuBinning::Range& r, Bool_t accEffCorrected) const
1013{
5376e016
CP		 1014 return TString::Format("MinvUS%s%s",r.AsString().Data(),
1015 accEffCorrected ? "_AccEffCorr" : "");
1016}
1017
16560e8e 1018
5376e016 1019//_____________________________________________________________________________
16560e8e 1020Double_t AliAnalysisMuMuMinv::GetAccxEff(Double_t pt,Double_t rapidity)
5376e016 1021{
16560e8e 1022 if (!fAccEffHisto)
1023 {
1024 AliError("ERROR: No AccxEff histo");
1025 return 0;
1026 }
1027 Int_t bin = fAccEffHisto->FindBin(pt,rapidity);
1028 Double_t accXeff = fAccEffHisto->GetBinContent(bin);
1029
1030 return accXeff;
1031}
1032
1033//_____________________________________________________________________________
1034Double_t AliAnalysisMuMuMinv::WeightDistribution(Double_t pt,Double_t rapidity)
1035{
1036 //Return a weight for a dimuon pt and y, which depend on the varied distributions.
1037 // FIXME: hard coded, find a clean way to fix the distribution parameters from outside
1038
1039 if (!HasMC() ) return 1.;
1040
1041 //================ p-Pb ==============//
1042 //value for input distribution: this is the nominal pt and y distribution
1043 Double_t paryPPB[2] = {1.0,0.174};
1044 Double_t parptPPB[3] = {1.0,0.0557,3.52};
1045
1046 Double_t paryHardPPB = 0.1344, parySoftPPB = 0.1971;
1047 Double_t par1ptHardPPB = 5.51e-2, par2ptHardPPB = 3.47,
1048 par1ptSoftPPB = 5.67e-2, par2ptSoftPPB = 3.68;
1049 //systematic 1: hardest in y x softest in pt
1050 Double_t pary1[2] = {1.0,paryHardPPB};
1051 Double_t parpt1[3] = {1.0,par1ptSoftPPB,par2ptSoftPPB};
1052 //systematic 2: hardest in y x hardest in pt
1053 Double_t pary2[2] = {1.0,paryHardPPB};
1054 Double_t parpt2[3] = {1.0,par1ptHardPPB,par2ptHardPPB};
1055 //systematic 3: softest in y x softest in pt
1056 Double_t pary3[2] = {1.0, parySoftPPB};
1057 Double_t parpt3[3] = {1.0,par1ptSoftPPB,par2ptSoftPPB};
1058 //systematic 4: softest in y x hardest in pt
1059 Double_t pary4[2] = {1.0,parySoftPPB};
1060 Double_t parpt4[3] = {1.0,par1ptHardPPB,par2ptHardPPB};
1061
1062 Double_t funcptvalPPB = powerLaw3Par(&pt,parptPPB);
1063 Double_t funcyvalPPB = normPol12Par(&rapidity,paryPPB);
1064
1065 //================ Pb-p ==============//
1066 //value for input distribution
1067 Double_t paryPBP[2] = {1.0,0.189};
1068 Double_t parptPBP[3] = {1.0,0.0592,3.92};
1069
1070 Double_t paryHardPBP = 0.1517, parySoftPBP = 0.2191;
1071 Double_t par1ptHardPBP = 5.58e-2, par2ptHardPBP = 3.83,
1072 par1ptSoftPBP = 5.59e-2, par2ptSoftPBP = 4.31;
1073 //systematic 5: hardest in y x softest in pt
1074 Double_t pary5[2] = {1.0,paryHardPBP};
1075 Double_t parpt5[3] = {1.0,par1ptSoftPBP,par2ptSoftPBP};
1076 //systematic 6: hardest in y x hardest in pt
1077 Double_t pary6[2] = {1.0,paryHardPBP};
1078 Double_t parpt6[3] = {1.0,par1ptHardPBP,par2ptHardPBP};
1079 //systematic 7: softest in y x softest in pt
1080 Double_t pary7[2] = {1.0, parySoftPBP};
1081 Double_t parpt7[3] = {1.0,par1ptSoftPBP,par2ptSoftPBP};
1082 //systematic 8: softest in y x hardest in pt
1083 Double_t pary8[2] = {1.0,parySoftPBP};
1084 Double_t parpt8[3] = {1.0,par1ptHardPBP,par2ptHardPBP};
1085
1086 Double_t funcptvalPBP = powerLaw3Par(&pt,parptPBP);
1087 Double_t funcyvalPBP = normPol12Par(&rapidity,paryPBP);
1088
1089 //================ pp ==============//
1090 //value for input distribution
1091 Double_t paryPP[2] = {3.0,0.514/3.};
1092 Double_t parptPP[3] = {1.0,0.0546,3.90};
1093
1094 Double_t paryHardPP = 0.4125/3., parySoftPP = 0.5958/3.;
1095 Double_t par1ptHardPP = 4.78e-2, par2ptHardPP = 3.65,//4.84e-2/3.45
1096 par1ptSoftPP = 6.12e-2, par2ptSoftPP = 4.31;//5.47e-2//4.29
1097 //systematic 9: hardest in y x softest in pt
1098 Double_t pary9[2] = {3.0,paryHardPP};
1099 Double_t parpt9[3] = {1.0,par1ptSoftPP,par2ptSoftPP};
1100 //systematic 10: hardest in y x hardest in pt
1101 Double_t pary10[2] = {3.0,paryHardPP};
1102 Double_t parpt10[3] = {1.0,par1ptHardPP,par2ptHardPP};
1103 //systematic 11: softest in y x softest in pt
1104 Double_t pary11[2] = {3.0, parySoftPP};
1105 Double_t parpt11[3] = {1.0,par1ptSoftPP,par2ptSoftPP};
1106 //systematic 12: softest in y x hardest in pt
1107 Double_t pary12[2] = {3.0,parySoftPP};
1108 Double_t parpt12[3] = {1.0,par1ptHardPP,par2ptHardPP};
1109
1110 Double_t funcptvalPP = powerLaw3Par(&pt,parptPP);
1111 Double_t funcyvalPP = normPol12Par(&rapidity,paryPP);
1112
1113 //______
1114 Double_t weight(1.),funcptsyst(0.),funcysyst(0.);
1115 switch ( fsystLevel )
1116 {
1117 case 0:
1118 weight = 1;
1119 break;
1120 case 1:
1121 funcptsyst = powerLaw3Par(&pt,parpt1);
1122 if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1123 else weight = 1;
1124 funcysyst = normPol12Par(&rapidity,pary1);
1125 if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1126 break;
1127 case 2:
1128 funcptsyst = powerLaw3Par(&pt,parpt2);
1129 if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1130 else weight = 1;
1131 funcysyst = normPol12Par(&rapidity,pary2);
1132 if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1133 break;
1134 case 3:
1135 funcptsyst = powerLaw3Par(&pt,parpt3);
1136 if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1137 else weight = 1;
1138 funcysyst = normPol12Par(&rapidity,pary3);
1139 if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1140 break;
1141 case 4:
1142 funcptsyst = powerLaw3Par(&pt,parpt4);
1143 if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1144 else weight = 1;
1145 funcysyst = normPol12Par(&rapidity,pary4);
1146 if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1147 break;
1148 case 5:
1149 funcptsyst = powerLaw3Par(&pt,parpt5);
1150 if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1151 else weight = 1;
1152 funcysyst = normPol12Par(&rapidity,pary5);
1153 if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1154 break;
1155 case 6:
1156 funcptsyst = powerLaw3Par(&pt,parpt6);
1157 if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1158 else weight = 1;
1159 funcysyst = normPol12Par(&rapidity,pary6);
1160 if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1161 break;
1162 case 7:
1163 funcptsyst = powerLaw3Par(&pt,parpt7);
1164 if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1165 else weight = 1;
1166 funcysyst = normPol12Par(&rapidity,pary7);
1167 if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1168 break;
1169 case 8:
1170 funcptsyst = powerLaw3Par(&pt,parpt8);
1171 if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1172 else weight = 1;
1173 funcysyst = normPol12Par(&rapidity,pary8);
1174 if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1175 break;
1176 case 9:
1177 funcptsyst = powerLaw3Par(&pt,parpt9);
1178 if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1179 else weight = 1;
1180 funcysyst = normPol12Par(&rapidity,pary9);
1181 if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1182 break;
1183 case 10:
1184 funcptsyst = powerLaw3Par(&pt,parpt10);
1185 if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1186 else weight = 1;
1187 funcysyst = normPol12Par(&rapidity,pary10);
1188 if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1189 break;
1190 case 11:
1191 funcptsyst = powerLaw3Par(&pt,parpt11);
1192 if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1193 else weight = 1;
1194 funcysyst = normPol12Par(&rapidity,pary11);
1195 if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1196 break;
1197 case 12:
1198 funcptsyst = powerLaw3Par(&pt,parpt12);
1199 if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1200 else weight = 1;
1201 funcysyst = normPol12Par(&rapidity,pary12);
1202 if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1203 break;
1204 }
1205
1206 return weight;
1207}
1208
1209//________________________________________________________________________
1210Double_t AliAnalysisMuMuMinv::powerLaw3Par(Double_t *x, Double_t *par)
1211{
1212 //3 parameters
1213 Double_t arg = 0;
1214
1215 arg = par[0]*x[0] / TMath::Power( 1 + par[1]*x[0]*x[0], par[2]);
1216
1217 return arg;
1218}
1219
1220
1221//________________________________________________________________________
1222Double_t AliAnalysisMuMuMinv::normPol12Par(Double_t *x, Double_t *par)
1223{
1224 //2 parameters
1225 Double_t arg1 = 0;
1226
1227 arg1 = par[0] * ( 1 + par[1]*x[0] );
1228
1229
1230 return arg1;
1231}
1232
1233//_____________________________________________________________________________
1234Bool_t AliAnalysisMuMuMinv::IsPtInRange(const AliVParticle& t1, const AliVParticle& t2, Double_t& ptmin, Double_t& ptmax) const
1235{
1236 /// Whether the pair passes the rapidity cut
5376e016
CP		 1237
1238 TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t1.P()*t1.P()));
1239 TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t2.P()*t2.P()));
1240
1241 TLorentzVector total(p1+p2);
1242
16560e8e 1243 Double_t pt = total.Pt();
1244
1245 return ( pt < ptmax && pt > ptmin );
1246}
5376e016 1247
16560e8e 1248//_____________________________________________________________________________
1249void AliAnalysisMuMuMinv::NameOfIsPtInRange(TString& name, Double_t& ptmin, Double_t& ptmax) const
1250{
1251 name.Form("PAIRPTIN%2.1f-%2.1f",ptmin,ptmax);
5376e016
CP		 1252}
1253
1254//_____________________________________________________________________________
16560e8e 1255Bool_t AliAnalysisMuMuMinv::IsRapidityInRange(const AliVParticle& t1, const AliVParticle& t2) const
5376e016 1256{
16560e8e 1257 TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t1.P()*t1.P()));
1258 TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t2.P()*t2.P()));
5376e016 1259
16560e8e 1260 TLorentzVector total(p1+p2);
1261
1262 Double_t y = total.Rapidity();
5376e016 1263
16560e8e 1264 return ( y < -2.5 && y > -4.0 );
1265}