]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muondep/AliAnalysisMuMu.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliAnalysisMuMu helper classes update
[u/mrichter/AliRoot.git] / PWG / muondep / AliAnalysisMuMu.cxx
CommitLineData
2f331ac9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17// $Id$
18
19#include "AliAnalysisMuMu.h"
cb690a12 20#include "AliAnalysisMuMuBase.h" //Just to have avaliable the MCInputPrefix() method
2f331ac9 21
a58729a5 22#include "AliAnalysisMuMuBinning.h"
cb690a12 23#include "AliAnalysisMuMuConfig.h"
81190958 24#include "AliAnalysisMuMuFnorm.h"
25#include "AliAnalysisMuMuGraphUtil.h"
26#include "AliAnalysisMuMuJpsiResult.h"
a58729a5 27#include "AliAnalysisMuMuSpectra.h"
2f331ac9 28#include "AliAnalysisTriggerScalers.h"
29#include "AliCounterCollection.h"
30#include "AliHistogramCollection.h"
31#include "AliLog.h"
a58729a5 32#include "AliMergeableCollection.h"
2f331ac9 33#include "Riostream.h"
34#include "TArrayL64.h"
a58729a5 35#include "TASImage.h"
36#include "TAxis.h"
37#include "TCanvas.h"
38#include "TColor.h"
39#include "TF1.h"
40#include "TFile.h"
41#include "TGraphErrors.h"
42#include "TGrid.h"
43#include "TH1.h"
44#include "TH2.h"
cb690a12 45#include "TProfile.h"
81190958 46#include "THashList.h"
a58729a5 47#include "TKey.h"
48#include "TLegend.h"
49#include "TLegendEntry.h"
50#include "TLine.h"
51#include "TList.h"
52#include "TMap.h"
53#include "TMath.h"
54#include "TObjArray.h"
55#include "TObjString.h"
56#include "TParameter.h"
57#include "TPaveText.h"
cb690a12 58#include "TPRegexp.h"
a58729a5 59#include "TROOT.h"
1afce1ce 60#include "TStopwatch.h"
a58729a5 61#include "TStyle.h"
62#include "TSystem.h"
63#include <cassert>
64#include <map>
65#include <set>
66#include <string>
cb690a12 67#include "TLatex.h"
2f331ac9 68
a58729a5 69ClassImp(AliAnalysisMuMu)
2f331ac9 70
2f331ac9 71//_____________________________________________________________________________
cb690a12 72AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName, const char* associatedSimFileName2, const char* beamYear) : TObject(),
2f331ac9 73fFilename(filename),
2f331ac9 74fCounterCollection(0x0),
a58729a5 75fBinning(0x0),
76fMergeableCollection(0x0),
77fRunNumbers(),
1afce1ce 78fCorrectionPerRun(0x0),
cb690a12 79fAssociatedSimulation(0x0),
80fAssociatedSimulation2(0x0),
81fParticleName(""),
82fConfig(0x0)
2f331ac9 83{
84 // ctor
85
a58729a5 86 GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
1afce1ce 87
cb690a12 88 if ( IsSimulation() )
1afce1ce 89 {
cb690a12 90 TString sFilename(filename);
91 if ( sFilename.Contains("JPSI",TString::kIgnoreCase) ) SetParticleName("JPsi"); //FIXME: DO a method for this
92 else if ( sFilename.Contains("PSIP",TString::kIgnoreCase) ) SetParticleName("PsiP");
93 else AliError("Unknown Particle Name in simulation");
94 }
1afce1ce 95
cb690a12 96 if ( fCounterCollection )
97 {
81190958 98 if ( strlen(associatedSimFileName) )
99 {
100 fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
cb690a12 101
102 TString sAssociatedSimFileName(associatedSimFileName);
103 if ( sAssociatedSimFileName.Contains("JPSI",TString::kIgnoreCase) ) fAssociatedSimulation->SetParticleName("JPsi");
104 else if ( sAssociatedSimFileName.Contains("PSIP",TString::kIgnoreCase) ) fAssociatedSimulation->SetParticleName("PsiP");
105 else AliError("Unknown Particle Name in associated simulation");
81190958 106 }
cb690a12 107
108 if ( strlen(associatedSimFileName2) )
109 {
110 fAssociatedSimulation2 = new AliAnalysisMuMu(associatedSimFileName2);
111
112 TString sAssociatedSimFileName2(associatedSimFileName2);
113 if ( sAssociatedSimFileName2.Contains("JPSI",TString::kIgnoreCase) ) fAssociatedSimulation2->SetParticleName("JPsi");
114 else if ( sAssociatedSimFileName2.Contains("PSIP",TString::kIgnoreCase) ) fAssociatedSimulation2->SetParticleName("PsiP");
115 else AliError("Unknown Particle Name in associated simulation 2");
116 }
117
118 fConfig = new AliAnalysisMuMuConfig(beamYear);
1afce1ce 119 }
2f331ac9 120}
121
122//_____________________________________________________________________________
123AliAnalysisMuMu::~AliAnalysisMuMu()
124{
125 // dtor
81190958 126
127 if ( fAssociatedSimulation )
128 {
129 fAssociatedSimulation->Update();
130 }
cb690a12 131 if ( fAssociatedSimulation2 )
132 {
133 fAssociatedSimulation2->Update();
134 }
81190958 135
136 Update();
137
2f331ac9 138 delete fCounterCollection;
a58729a5 139 delete fBinning;
140 delete fMergeableCollection;
141 delete fCorrectionPerRun;
1afce1ce 142 delete fAssociatedSimulation;
cb690a12 143 delete fAssociatedSimulation2;
144 delete fConfig;
2f331ac9 145}
146
147//_____________________________________________________________________________
2754fd07 148void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
149 ULong64_t* totalNmb,
150 ULong64_t* totalNmsl,
151 ULong64_t* totalNmul)
2f331ac9 152{
153 // Report of some basic numbers, like number of MB and MUON triggers,
154 // both before and after physics selection, and comparison with
155 // the total number of such triggers (taken from the OCDB scalers)
156 // if requested.
157 //
158 // filename is assumed to be a root filecontaining a list containing
159 // an AliCounterCollection (or directly an AliCounterCollection)
160 //
161 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
162 //
163
a58729a5 164 if (!fMergeableCollection || !fCounterCollection) return;
2f331ac9 165
166 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
167 TIter nextRun(runs);
168
169 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
170 TIter nextTrigger(triggers);
171
172 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
173
2754fd07 174 Bool_t doPS = (events->FindObject("PSALL") != 0x0);
175
2f331ac9 176 TObjString* srun;
177 TObjString* strigger;
178
2754fd07 179 ULong64_t localNmb(0);
180 ULong64_t localNmsl(0);
181 ULong64_t localNmul(0);
182
183 if ( totalNmb) *totalNmb = 0;
184 if ( totalNmsl) *totalNmsl = 0;
185 if ( totalNmul ) *totalNmul = 0;
2f331ac9 186
187 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
188 {
189 std::cout << Form("RUN %09d ",srun->String().Atoi());
190
191 TString details;
192 ULong64_t nmb(0);
193 ULong64_t nmsl(0);
194 ULong64_t nmul(0);
195
196 nextTrigger.Reset();
197
2754fd07 198 Int_t nofPS(0);
199
2f331ac9 200 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
201 {
2754fd07 202
cb690a12 203 if ( !Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,IsSimulation()).Contains(strigger->String().Data()) &&
204 !Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation()).Contains(strigger->String().Data()) &&
205 !Config()->GetList(AliAnalysisMuMuConfig::kMuonTriggerList,IsSimulation()).Contains(strigger->String().Data()) ) continue;
2754fd07 206
81dfe194 207 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
208 strigger->String().Data(),"ALL",srun->String().Atoi())));
2f331ac9 209
2754fd07 210 details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
2f331ac9 211
212
81dfe194 213 ULong64_t nps = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
214 strigger->String().Data(),"PSALL",srun->String().Atoi())));
2754fd07 215
216 if ( doPS )
2f331ac9 217 {
2754fd07 218 details += TString::Format(" PS %5.1f %%",nps*100.0/n);
2f331ac9 219 }
220
2754fd07 221 if (nps)
222 {
223 ++nofPS;
224 }
225
cb690a12 226 if ( Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,IsSimulation()).Contains(strigger->String()) )
2f331ac9 227 {
228 nmb += n;
2754fd07 229 if ( totalNmb) (*totalNmb) += n;
230 localNmb += n;
2f331ac9 231 }
cb690a12 232 else if ( Config()->GetList(AliAnalysisMuMuConfig::kMuonTriggerList,IsSimulation()).Contains(strigger->String()) )
2f331ac9 233 {
234 nmsl += n;
2754fd07 235 if ( totalNmsl) (*totalNmsl) += n;
236 localNmsl += n;
2f331ac9 237 }
cb690a12 238 else if ( Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation()).Contains(strigger->String()) )
2f331ac9 239 {
240 nmul += n;
2754fd07 241 if ( totalNmul ) (*totalNmul) += n;
242 localNmul += n;
2f331ac9 243 }
244 }
245
2754fd07 246 std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
247 nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
2f331ac9 248
249 if ( detailTriggers )
250 {
251 std::cout << details.Data();
252 }
253 std::cout << std::endl;
254 }
255
2754fd07 256 if ( !totalNmul && !totalNmsl && !totalNmb )
257 {
258 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
259 localNmb,localNmsl,localNmul) << std::endl;
260 }
2f331ac9 261
2754fd07 262 delete runs;
2f331ac9 263 delete triggers;
264 delete events;
265}
266
2754fd07 267
a58729a5 268
1afce1ce 269//_____________________________________________________________________________
270void AliAnalysisMuMu::CleanAllSpectra()
271{
272 /// Delete all the spectra we may have
273
cb690a12 274 OC()->RemoveByType("AliAnalysisMuMuSpectra");
1afce1ce 275 Update();
276}
a58729a5 277
a58729a5 278
2f331ac9 279//_____________________________________________________________________________
280TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
281 const char* backgroundresults)
282{
283 TFile* fjpsi = FileOpen(jpsiresults);
284 TFile* fbck = FileOpen(backgroundresults);
285
286 if (!fjpsi || !fbck) return 0x0;
287
288 TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
289
290 std::vector<std::string> checks;
291
292 checks.push_back("muminus-CMUU7-B-NOPF-ALLNOTRD");
293 checks.push_back("muplus-CMUU7-B-NOPF-ALLNOTRD");
294 checks.push_back("muminus-CMUSH7-B-NOPF-MUON");
295 checks.push_back("muplus-CMUSH7-B-NOPF-MUON");
296
297 if (!gjpsi) return 0x0;
298
299 TObjArray* a = new TObjArray;
300 a->SetOwner(kTRUE);
301
302 for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
303 {
304
305 TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
306
307 if (!gback) continue;
308
309 if ( gjpsi->GetN() != gback->GetN() )
310 {
311 AliErrorClass("graphs have different number of points !");
312 continue;
313 }
314
315 TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
316
317 for ( int i = 0; i < gjpsi->GetN(); ++i )
318 {
319 double r1,r2,y1,y2;
320
321 gjpsi->GetPoint(i,r1,y1);
322 gback->GetPoint(i,r2,y2);
323
324 if ( r1 != r2 )
325 {
326 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
327 continue;
328 }
329
330 g->SetPoint(i,y2,y1);
331 // g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
332 }
333
334 g->SetMarkerStyle(25+j);
335 g->SetMarkerSize(1.2);
336 if (j==0)
337 {
338 g->Draw("ap");
339 }
340 else
341 {
342 g->Draw("p");
343 }
344 g->SetLineColor(j+1);
345 g->SetMarkerColor(j+1);
346 g->SetName(checks[j].c_str());
347 a->AddLast(g);
348 }
349
350 return a;
351}
352
353//_____________________________________________________________________________
354TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
355 const char* simjpsiresults)
356{
357 TFile* freal = FileOpen(realjpsiresults);
358 TFile* fsim = FileOpen(simjpsiresults);
359
360 if (!freal || !fsim) return 0x0;
361
362 TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
363 TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
364
365 TObjArray* a = new TObjArray;
366 a->SetOwner(kTRUE);
367
368 if ( greal->GetN() != gsim->GetN() )
369 {
370 AliErrorClass("graphs have different number of points !");
371 return 0x0;
372 }
373
374 TGraphErrors* g = new TGraphErrors(greal->GetN());
375 TGraphErrors* gratio = new TGraphErrors(greal->GetN());
376
377 for ( int i = 0; i < greal->GetN(); ++i )
378 {
379 double r1,r2,y1,y2;
380
381 greal->GetPoint(i,r1,y1);
382 gsim->GetPoint(i,r2,y2);
383
384 if ( r1 != r2 )
385 {
386 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
387 continue;
388 }
389
390 double ratio(0.0);
391
392 if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
393 {
394 g->SetPoint(i,0,0);
395 g->SetPointError(i,0,0);
396 }
397 else
398 {
399 g->SetPoint(i,y2,y1);
400 g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
401 ratio = y2/y1;
402 }
403 gratio->SetPoint(i,r1,ratio);
404 }
405
406 g->SetMarkerStyle(25);
407 g->SetMarkerSize(1.2);
408
409 new TCanvas;
410
411 g->Draw("ap");
412
413 g->SetLineColor(1);
414 g->SetMarkerColor(1);
415 g->SetName("jpsipercmuurealvssim");
416
417 new TCanvas;
418
419 greal->Draw("alp");
420 gsim->SetLineColor(4);
421
422 gsim->Draw("lp");
423
424 new TCanvas;
425 gratio->Draw("alp");
426
427 return g;
428}
429
430//_____________________________________________________________________________
cb690a12 431AliAnalysisMuMuConfig* AliAnalysisMuMu::Config()
432{
433 /// Return the configuration
434 return fConfig;
435}
436
437//_____________________________________________________________________________
438void AliAnalysisMuMu::DrawMinv(const char* type,
439 const char* particle,
440 const char* trigger,
441 const char* eventType,
442 const char* pairCut,
443 const char* centrality,
444 const char* subresultname,
445 const char* flavour) const
2f331ac9 446{
cb690a12 447 /// Draw minv spectra for binning of given type
448
449 if (!OC() || !BIN()) return;
2f331ac9 450
cb690a12 451 TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
452 if (!bins)
453 {
454 AliError(Form("Could not get %s bins",type));
455 return;
456 }
2f331ac9 457
cb690a12 458 Double_t xmin(-1);
459 Double_t xmax(-1);
2f331ac9 460
cb690a12 461 TString sparticle(particle);
462 if ( sparticle=="PSI" )
463 {
464 xmin = 2;
465 xmax = 6;
466 }
2f331ac9 467
cb690a12 468 Int_t nx(1);
469 Int_t ny(1);
2f331ac9 470
cb690a12 471 Int_t n = bins->GetEntries();
2f331ac9 472
cb690a12 473 if ( n == 2 )
474 {
475 nx = 2;
476 }
477 else if ( n > 2 )
478 {
479 ny = TMath::Nint(TMath::Sqrt(n));
480 nx = n/ny;
481 }
2f331ac9 482
cb690a12 483 TString stype(type);
484 stype.ToUpper();
485
486 TString spectraName(Form("/%s/%s/%s/%s/%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data()));
487
488 if ( strlen(flavour))
489 {
490 spectraName += "-";
491 spectraName += flavour;
492 }
493
494 AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(spectraName.Data()));
495
496 AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
2f331ac9 497
cb690a12 498 TObjArray* spectraBins(0x0);
499 if ( spectra )
500 {
501 spectraBins = spectra->BinContentArray();
502 }
2f331ac9 503
cb690a12 504 TCanvas* c = new TCanvas;
505 c->Divide(nx,ny);
506 c->Draw();
507 gStyle->SetOptFit(1112);
508
509 c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
510 (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
511
512
513 TIter next(bins);
514 AliAnalysisMuMuBinning::Range* r;
515 Int_t ci(0);
516
517 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
2f331ac9 518 {
cb690a12 519 TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
520
521 AliDebug(1,name.Data());
2f331ac9 522
cb690a12 523 AliAnalysisMuMuJpsiResult* spectraBin(0x0);
2f331ac9 524
cb690a12 525 if ( spectraBins )
2f331ac9 526 {
cb690a12 527 AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
a58729a5 528
cb690a12 529 spectraBin = static_cast<AliAnalysisMuMuJpsiResult*>(sr->SubResult(subresultname));
a58729a5 530
cb690a12 531 AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
2f331ac9 532 }
533
cb690a12 534 TH1* h = OC()->Histo(name.Data());
535
536 if ( spectraBin )
2f331ac9 537 {
cb690a12 538 h = spectraBin->Histo();
2f331ac9 539 }
2f331ac9 540
cb690a12 541 if (h)
2f331ac9 542 {
cb690a12 543 ++ci;
544 c->cd(ci);
545 gPad->SetLogy();
546 if (xmin>0)
2f331ac9 547 {
cb690a12 548 h->GetXaxis()->SetRangeUser(xmin,xmax);
2f331ac9 549 }
cb690a12 550 h->Draw("histes");
2f331ac9 551
cb690a12 552 TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
553 if (f1)
2f331ac9 554 {
cb690a12 555 f1->Draw("same");
2f331ac9 556 }
557
cb690a12 558 gPad->Modified();
559 gPad->Update();
2f331ac9 560
cb690a12 561 TObject* stats = h->FindObject("stats");
562 if (stats)
2f331ac9 563 {
cb690a12 564 stats->Draw("same");
2f331ac9 565 }
566 }
2f331ac9 567 }
568
cb690a12 569 delete bins;
570}
571
572//_____________________________________________________________________________
573void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
574{
575 /// Draw minv spectra for binning of given type
576
577// AliWarning("Reimplement me!");
2f331ac9 578
cb690a12 579 if (!fConfig)
2f331ac9 580 {
cb690a12 581 AliError("No configuration available yet. Don't know what to draw");
582 return;
2f331ac9 583 }
584
cb690a12 585 DrawMinv(type,particle,
586 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation())),
587 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation())),
588 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation())),
589 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation())),
590 subresultname,
591 flavour);
592// First(Config()->kEventSelectionList()).Data(),
593// First(Config()->kPairSelectionList()).Data(),
594// First(Config()->kCentralitySelectionList()).Data(),
595// subresultname,
596// flavour);
2f331ac9 597}
598
cb690a12 599//___________________________________________________________________
600void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
2754fd07 601{
cb690a12 602 // Actions in reponse to mouse button events.
2754fd07 603
cb690a12 604 TCanvas* c = static_cast<TCanvas*>(gTQSender);
605 TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
606 if (!pad) return;
2754fd07 607
cb690a12 608// if ((event == kButton1Down) ||
609 if (event == kButton1Double)
2754fd07 610 {
2754fd07 611
cb690a12 612// Float_t x = pad->AbsPixeltoX(px);
613// Float_t y = pad->AbsPixeltoY(py);
614// x = pad->PadtoX(x);
615// y = pad->PadtoY(y);
616
617// std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
2754fd07 618
cb690a12 619 if ( sel && sel->InheritsFrom("TH1") )
620 {
621 TCanvas* clocal = new TCanvas;
622 clocal->SetLogy();
623 clocal->Draw();
624 sel->Draw();
625 }
626 else
2754fd07 627 {
cb690a12 628 TList* list = pad->GetListOfPrimitives();
629 TIter next(list);
630 TObject* h;
2754fd07 631
cb690a12 632 while ( ( h = next() ) )
633 {
634 if ( h->InheritsFrom("TH1") )
635 {
636 TCanvas* clocal = new TCanvas;
637 clocal->SetLogy();
638 clocal->Draw();
639 h->Draw();
640 break;
641 }
642 }
2754fd07 643
2754fd07 644 }
2754fd07 645
cb690a12 646// std::cout << std::endl;
a58729a5 647
cb690a12 648 pad->Modified();
649 }
2754fd07 650
cb690a12 651}
a58729a5 652
cb690a12 653//_____________________________________________________________________________
654TString
655AliAnalysisMuMu::ExpandPathName(const char* file)
656{
657 // An expand method that lives alien URL as they are
658 TString sfile;
2754fd07 659
cb690a12 660 if ( !sfile.BeginsWith("alien://") )
661 {
662 return gSystem->ExpandPathName(file);
663 }
664 else
665 {
666 if (!gGrid) TGrid::Connect("alien://");
667 if (!gGrid) return "";
668 }
2754fd07 669
cb690a12 670 return file;
2754fd07 671}
672
2f331ac9 673//_____________________________________________________________________________
cb690a12 674void AliAnalysisMuMu::TwikiOutputFnorm(const char* series) const
2f331ac9 675{
cb690a12 676 // make a twiki-compatible output of the Fnorm factor(s)
677 TObjArray* what = TString(series).Tokenize(",");
678 TObjString* s;
679 TObjArray graphs;
680 TIter next(what);
681
682 std::cout << "| *Run* |";
683 while ( ( s = static_cast<TObjString*>(next())) )
684 {
685 TGraph* g = static_cast<TGraph*>(OC()->GetObject(Form("/FNORM/GRAPHS/%s",s->String().Data())));
686 if (!g)
687 {
688 AliError(Form("Could not find graph for %s",s->String().Data()));
689 continue;
690 }
691 std::cout << " *" << s->String().Data();
692 if ( s->String().BeginsWith("RelDif") ) std::cout << " %";
693 std::cout << "*|";
694 graphs.Add(g);
695 }
81190958 696
cb690a12 697 std::cout << std::endl;
2f331ac9 698
cb690a12 699 TGraphErrors* g0 = static_cast<TGraphErrors*>(graphs.First());
700 if (!g0) return;
2f331ac9 701
cb690a12 702 for ( Int_t i = 0; i < g0->GetN(); ++i )
2f331ac9 703 {
cb690a12 704 TString msg;
81190958 705
cb690a12 706 msg.Form("|%6d|",TMath::Nint(g0->GetX()[i]));
81190958 707
cb690a12 708 for ( Int_t j = 0; j < graphs.GetEntries(); ++j )
81190958 709 {
cb690a12 710 TGraphErrors* g = static_cast<TGraphErrors*>(graphs.At(j));
711
712 msg += TString::Format(" %6.2f +- %6.2f |",g->GetY()[i],g->GetEY()[i]);
81190958 713 }
81190958 714
cb690a12 715 std::cout << msg.Data() << std::endl;
2f331ac9 716 }
717
cb690a12 718 next.Reset();
719
720 std::cout << "|*Weigthed mean (*)*|";
721
722 AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
2f331ac9 723
cb690a12 724 if (!r)
2f331ac9 725 {
cb690a12 726 AliError("Could not find Fnorm result !");
727 return;
728 }
729
730
731 while ( ( s = static_cast<TObjString*>(next())) )
732 {
733 TString var("Fnorm");
734 TString unit;
735
736 if ( s->String().BeginsWith("Fnorm") )
2f331ac9 737 {
cb690a12 738 r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
2f331ac9 739 }
cb690a12 740 else if ( s->String().BeginsWith("RelDif") )
2f331ac9 741 {
cb690a12 742 r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/RelDif"));
743 unit = "%";
2f331ac9 744 }
cb690a12 745
746 r->Exclude("*");
747 r->Include(s->String().Data());
748
749 std::cout << Form(" * %5.2f +- %5.2f %s * |",
750 r->GetValue(var.Data()),
751 r->GetErrorStat(var.Data()),
752 unit.Data());
753 }
754
755 next.Reset();
756
757 std::cout << std::endl;
758
759 std::cout << "|*RMS*|";
760
761 while ( ( s = static_cast<TObjString*>(next())) )
762 {
763 TString var("Fnorm");
764
765 if ( s->String().BeginsWith("Fnorm") )
2f331ac9 766 {
cb690a12 767 r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
2f331ac9 768 }
cb690a12 769 else if ( s->String().BeginsWith("RelDif") )
2f331ac9 770 {
cb690a12 771 r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/RelDif"));
2f331ac9 772 }
cb690a12 773
774 r->Exclude("*");
775 r->Include(s->String().Data());
776
777 Double_t d = 100.0*r->GetRMS(var.Data())/r->GetValue(var.Data());
778
779 std::cout << Form(" * %5.2f (%5.2f %%) * |",
780 r->GetRMS(var.Data()),d);
2f331ac9 781 }
782
cb690a12 783 std::cout << std::endl;
784 std::cout << "(*) weight is the number of CMUL7-B-NOPF-MUON triggers (physics-selected and pile-up corrected) in each run" << std::endl;
2f331ac9 785
cb690a12 786 delete what;
2f331ac9 787}
788
a58729a5 789//_____________________________________________________________________________
cb690a12 790TFile*
791AliAnalysisMuMu::FileOpen(const char* file)
a58729a5 792{
cb690a12 793 // Open a file after expansion of its name
a58729a5 794
cb690a12 795 return TFile::Open(ExpandPathName(file).Data());
796}
797
798//_____________________________________________________________________________
799TString AliAnalysisMuMu::First(const TString& list) const
800{
801 TObjArray* a = list.Tokenize(",");
802 if ( a->GetLast() < 0 ) return "";
1afce1ce 803
cb690a12 804 TString rv = static_cast<TObjString*>(a->First())->String();
1afce1ce 805
cb690a12 806 delete a;
1afce1ce 807
cb690a12 808 return rv;
809}
810
811//_____________________________________________________________________________
812AliAnalysisMuMuSpectra*
813AliAnalysisMuMu::FitParticle(const char* particle,
814 const char* trigger,
815 const char* eventType,
816 const char* pairCut,
817 const char* centrality,
818 const AliAnalysisMuMuBinning& binning,
819 const char* spectraType,
820 Bool_t corrected)
821{
822 // Fit the minv/mpt spectra to find the given particle
823 // Returns an array of AliAnalysisMuMuResult objects
1afce1ce 824
cb690a12 825 TProfile::Approximate(); //To avoid bins with error=0 due to low statstics
826
827 static int n(0);
828
829 TObjArray* bins = binning.CreateBinObjArray(particle);
830 if (!bins)
1afce1ce 831 {
cb690a12 832 AliError(Form("Did not get any bin for particle %s",particle));
833 return 0x0;
1afce1ce 834 }
cb690a12 835
836 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
837 if ( !triggers->FindObject(trigger) )
1afce1ce 838 {
cb690a12 839 AliError(Form("Did not find trigger %s",trigger));
840 delete bins;
841 delete triggers;
842 return 0x0;
1afce1ce 843 }
cb690a12 844 delete triggers;
1afce1ce 845
cb690a12 846 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
847 if ( !events->FindObject(eventType) )
81190958 848 {
cb690a12 849 AliError(Form("Did not find eventType %s",eventType));
850 delete bins;
851 delete events;
852 return 0x0;
81190958 853 }
cb690a12 854 delete events;
1afce1ce 855
cb690a12 856 Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
857
858 if (ntrigger<=0)
1afce1ce 859 {
cb690a12 860 AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
861 delete bins;
862 return 0x0;
1afce1ce 863 }
cb690a12 864
865 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
866 Int_t nruns = runs->GetEntries();
867 delete runs;
1afce1ce 868
1afce1ce 869
cb690a12 870 TString id(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut));
871
872// binning.Print();
873
874 AliAnalysisMuMuSpectra* spectra(0x0);
1afce1ce 875
cb690a12 876 AliAnalysisMuMuBinning::Range* bin;
a58729a5 877 TIter next(bins);
a58729a5 878
cb690a12 879 TObjArray* fitTypeArray = Config()->GetList(AliAnalysisMuMuConfig::kFitTypeList,IsSimulation()).Tokenize(",");
880 TIter nextFitType(fitTypeArray);
881 TObjString* fitType;
882 TString flavour;
883 TString sSpectraType(spectraType);
884
885 TString spectraName(binning.GetName());
886 if ( flavour.Length() > 0 )
a58729a5 887 {
cb690a12 888 spectraName += "-";
889 spectraName += flavour;
890 }
891 if ( corrected )
892 {
893 spectraName += "-";
894 spectraName += "AccEffCorr";
895 }
896
897// Int_t binN(0);
898 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
899 {
900 TString hname;
901 if (!sSpectraType.CompareTo("minv"))
1afce1ce 902 {
cb690a12 903 hname = corrected ? Form("MinvUS%s_AccEffCorr",bin->AsString().Data()) : Form("MinvUS%s",bin->AsString().Data());
904 }
905 else if (!sSpectraType.CompareTo("mpt"))
906 {
907 hname = corrected ? Form("MeanPtVsMinvUS%s_AccEffCorr",bin->AsString().Data()) : Form("MeanPtVsMinvUS%s",bin->AsString().Data());
908 }
909 else
910 {
911 AliError("Wrong spectra type choice: Posibilities are: 'minv' or 'mpt' ");
912 return 0x0;
1afce1ce 913 }
914
cb690a12 915 TString isCorr(corrected ? " AccEffCorr " : " ");
916 std::cout << "---------------------------------//---------------------------------" << std::endl;
917 std::cout << "Fitting" << isCorr.Data() << sSpectraType.Data() << " spectra in " << id.Data() << std::endl;
1afce1ce 918
cb690a12 919 TH1* histo = OC()->Histo(id.Data(),hname.Data());
920
921 if (!histo)
1afce1ce 922 {
cb690a12 923// if (!fBinning && bin->IsIntegrated() )
924// {
925// // old file, we only had MinvUSPt
926// hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
927// }
928//
929// if (!hminv)
930// {
931 AliError(Form("Could not find histo %s",hname.Data()));
932 continue;
933// }
1afce1ce 934 }
935
cb690a12 936 histo = static_cast<TH1*>(histo->Clone(Form("%s%d",sSpectraType.Data(),n++)));
937
938 const char* particleTmp = IsSimulation() ? GetParticleName() : "JPsi"; //At some point particleTmp should become particle (but for now particle is always = "psi")
939
940 TString sparticleTmp(particleTmp);
941
942 AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(particleTmp,
943 *histo, // Result for current bin
944 trigger,
945 eventType,
946 pairCut,
947 centrality,
948 *bin);
949
950 r->SetNofTriggers(ntrigger);
951 r->SetNofRuns(nruns);
952
953 nextFitType.Reset();
954
955 Int_t added(0);
956
957 while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
a58729a5 958 {
cb690a12 959 // In this loop we create a Subresult for each fit inside the Result for current bin (AddFit will do)
960
961 TString sFitType(fitType->String());
962
963 if ( !sFitType.Contains(sSpectraType.Data()) ) continue;
964
965 AliDebug(1,Form("<<<<<< fitType=%s bin=%s",sFitType.Data(),bin->Flavour().Data()));
966
967 std::cout << "" << std::endl;
968 std::cout << "---------------" << "Fit " << added + 1 << "------------------" << std::endl;
969 std::cout << "Fitting " << hname.Data() << " with " << sFitType.Data() << std::endl;
970 std::cout << "" << std::endl;
971
972 if ( sFitType.Contains("mctails",TString::kIgnoreCase) ) //FIXME: Find a univoque way to determine the correctly the fit type
1afce1ce 973 {
cb690a12 974 TString sbin = bin->AsString();
975 TString spectraMCName = spectraName;
976 AliAnalysisMuMuBinning::Range* binMC = bin;
977
978 if ((sbin.Contains("MULT") || sbin.Contains("NCH") || sbin.Contains("DNCHDETA") || sbin.Contains("V0A") || sbin.Contains("V0ACENT") || sbin.Contains("NTRCORR")|| sbin.Contains("RELNTRCORR")) && !sbin.Contains("NTRCORRPT") && !sbin.Contains("NTRCORRY"))
979 {
980 //-------has to have a better way to do it
981 AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
982 b->AddBin("psi","INTEGRATED");
983
984 binMC = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
985
986 spectraMCName = b->GetName();
987 delete b;
988
989// if ( flavour.Length() > 0 ) //Commented cause we set no flavour in "INTEGRATED" bin in the analysis task
990// {
991// spectraMCName += "-";
992// spectraMCName += flavour;
993// }
994 if ( corrected )
995 {
996 spectraMCName += "-";
997 spectraMCName += "AccEffCorr";
998 }
999 //-----------
1000
1001 }
1002// if( sbin.Contains("NTRCORRPT") || !sbin.Contains("NTRCORRY") )
1003// {
1004// spectraMCName.Remove(4,7);
1005//
1006// if ( spectraMCName.Contains("PT") )
1007// {
1008// AliAnalysisMuMuBinning* b = SIM()->BIN();
1009//
1010// TObjArray* binsPt = b->CreateBinObjArray(particle,"PT","");
1011//
1012// Int_t nEntrBinMC = binsPt->GetEntries();
1013//
1014//// binsPt->Print();
1015//
1016// binMC = static_cast<AliAnalysisMuMuBinning::Range*>(binsPt->At(binN));
1017//
1018// if ( binN == nEntrBinMC - 1 ) binN = -1;
1019//
1020// }
1021//
1022// }
1023
1024 //par = GetCB2Tails(*binInt,"MCTAILS",eventType,trigger,pairCut,corrected);Why I was taking the tails from the fitted data spectra?
1025 GetParametersFromMC(sFitType,Form("/%s/%s",centrality,pairCut),spectraMCName.Data(),binMC);
1026
1027 if (sFitType.Length()>0)
1028 {
1029 added += ( r->AddFit(sFitType.Data()) == kTRUE );
1030 }
1afce1ce 1031 }
1afce1ce 1032
cb690a12 1033 else if ( sFitType.Contains("mpt",TString::kIgnoreCase) && !sFitType.Contains("minv",TString::kIgnoreCase) )
1afce1ce 1034 {
cb690a12 1035 std::cout << "++The Minv parameters will be taken from " << spectraName.Data() << std::endl;
1036 std::cout << "" << std::endl;
1037
1038 AliAnalysisMuMuSpectra* minvSpectra = dynamic_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(id.Data(),spectraName.Data()));
1039
1040 if ( !minvSpectra )
1041 {
1042 AliError(Form("Cannot fit mean pt: could not get the minv spectra for %s",id.Data()));
1043 continue;//return 0x0;
1044 }
1045
1046 AliAnalysisMuMuJpsiResult* minvResult = static_cast<AliAnalysisMuMuJpsiResult*>(minvSpectra->GetResultForBin(*bin));
1047
1048 if ( !minvResult )
1049 {
1050 AliError(Form("Cannot fit mean pt: could not get the minv result for bin %s in %s",bin->AsString().Data(),id.Data()));
1051 continue; //return 0x0;
1052 }
1053
1054 TObjArray* minvSubResults = minvResult->SubResults();
1055 TIter nextSubResult(minvSubResults);
1056 AliAnalysisMuMuJpsiResult* fitMinv;
1057 TString subResultName;
1058
1059 Int_t nSubFit(0);
1060 while ( ( fitMinv = static_cast<AliAnalysisMuMuJpsiResult*>(nextSubResult())) )
1061 {
1062 std::cout << "" << std::endl;
1063 std::cout << " /-- SubFit " << nSubFit + 1 << " --/ " << std::endl;
1064 std::cout << "" << std::endl;
1065
1066 TString fitMinvName(fitMinv->GetName());
1067 fitMinvName.Remove(fitMinvName.First("_"),fitMinvName.Sizeof()-fitMinvName.First("_"));
1068
1069// std::cout << fitMinvName.Data() << " ; " << fitMinv->GetName() << std::endl;
1070// std::cout << "" << std::endl;
1071
1072 if ( !sFitType.Contains(fitMinvName) ) continue; //FIXME: Ambiguous, i.e. NA60NEWPOL2EXP & NA60NEWPOL2 (now its ok cause only VWG and POL2EXP are used, but care)
1073
1074 TString sMinvFitType(sFitType);
1075
1076 GetParametersFromResult(sMinvFitType,fitMinv);//FIXME: Think about if this is necessary
1077
1078 added += ( r->AddFit(sMinvFitType.Data()) == kTRUE );
1079
1080 nSubFit++;
1081 }
1afce1ce 1082 }
1083
cb690a12 1084 else if ( sFitType.Contains("minv&mpt",TString::kIgnoreCase) )
1085 {
1086 AliWarning("Implement here the method to do the combined minv mpt fits");
1087 //FIXME: Shall we use the fitType or spectraType to choose to perform combined fits? Cause we have to check what kind of object is returned by the combined fit in order to decide if we put it in a different spectra(spectraType would be the flag,and here we should update the spectraName) or as a subresult(fitType in this case)
1088 }
1afce1ce 1089
cb690a12 1090 else
1afce1ce 1091 {
cb690a12 1092 if ( sFitType.Contains("PSICB2",TString::kIgnoreCase) || sFitType.Contains("PSINA60NEW",TString::kIgnoreCase)) std::cout << "+Free tails fit... " << std::endl;
1093 else if ( sFitType.Contains("PSICOUNT",TString::kIgnoreCase) ) std::cout << Form("+Just counting %s...",GetParticleName()) << std::endl;
1094 else std::cout << "+Using predefined tails... " << std::endl;
1095
1096 if ( sFitType.Contains("minvJPsi") && !sparticleTmp.Contains("JPsi") )
1097 {
1098 std::cout << "This fitting funtion is set to fit JPsi: Skipping fit..." << std::endl;
1099 continue;
1100 }
1101 if ( sFitType.Contains("minvPsiP") && !sparticleTmp.Contains("PsiP") )
1102 {
1103 std::cout << "This fitting funtion is set to fit PsiP: Skipping fit..." << std::endl;
1104 continue;
1105 }
1106
1107 added += ( r->AddFit(sFitType.Data()) == kTRUE );
1afce1ce 1108 }
cb690a12 1109
1110 std::cout << "-------------------------------------" << std::endl;
1111 std::cout << "" << std::endl;
a58729a5 1112 }
1afce1ce 1113
cb690a12 1114 if ( !added ) continue;
1afce1ce 1115
cb690a12 1116// TObjArray* a = r->SubResults(); // TEST
1117// a->Print();
1afce1ce 1118
cb690a12 1119 flavour = bin->Flavour();
1120
1121 if (!spectra)
1afce1ce 1122 {
cb690a12 1123 TString spectraSaveName = spectraName;
1124 if ( !sSpectraType.CompareTo("mpt") )
1afce1ce 1125 {
cb690a12 1126 spectraSaveName += "-";
1127 spectraSaveName += "MeanPtVsMinvUS";
1afce1ce 1128 }
1129
cb690a12 1130 spectra = new AliAnalysisMuMuSpectra(spectraSaveName.Data());
1afce1ce 1131 }
cb690a12 1132
1133 Bool_t adoptOk = spectra->AdoptResult(*bin,r); // We adopt the Result for current bin into the spectra
1134
1135 if ( adoptOk ) std::cout << "Result " << r->GetName() << " adopted in spectra " << spectra->GetName() << std::endl;
1136 else AliError(Form("Error adopting result %s in spectra %s",r->GetName(),spectra->GetName()));
1137
1138
1139 if ( IsSimulation() )
1140 {
1141 SetNofInputParticles(*r);
1142 }
1143
1144// binN++;
1145 } // loop on bins
1146
1147 delete fitTypeArray;
1148 delete bins;
1149
1150 return spectra;
1151}
1152
1153
1154////_____________________________________________________________________________
1155//AliAnalysisMuMuSpectra*
1156//AliAnalysisMuMu::FitMeanPt(const char* particle,
1157// const char* trigger,
1158// const char* eventType,
1159// const char* pairCut,
1160// const char* centrality,
1161// const AliAnalysisMuMuBinning& binning,const AliAnalysisMuMuSpectra& spectraMinv,Bool_t corrected)
1162//{
1163// // Fit the dimuon mean pt spectra to find the particle mean pt
1164// // Returns an array of AliAnalysisMuMuResult objects
1165//
1166// TObjArray* bins = binning.CreateBinObjArray(particle);
1167// if (!bins)
1168// {
1169// AliError(Form("Did not get any bin for particle %s",particle));
1170// return 0x0;
1171// }
1172//
1173// // AliAnalysisMuMuBinning* binningMinv(0x0);
1174// // TObjArray* binsMinv(0x0);
1175// // AliAnalysisMuMuBinning::Range* binInt(0x0);
1176// // AliAnalysisMuMuBinning::Range* testbin = static_cast<AliAnalysisMuMuBinning::Range*>(bins->At(0));
1177// // TString sbins = testbin->AsString();
1178//
1179// // if ( sbins.Contains("MULT")) //For the multiplicity bins we will use the tails from the integrated Minv spectra
1180// // {
1181// // binningMinv = spectraMinv.Binning();
1182// // binsMinv = binningMinv->CreateBinObjArray(particle);
1183// // binInt = static_cast<AliAnalysisMuMuBinning::Range*>(binsMinv->At(0));
1184// // }
1185//
1186//
1187// TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1188// if ( !triggers->FindObject(trigger) )
1189// {
1190// AliDebug(1,Form("Did not find trigger %s",trigger));
1191// delete bins;
1192// delete triggers;
1193// return 0x0;
1194// }
1195// delete triggers;
1196//
1197// TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1198// if ( !events->FindObject(eventType) )
1199// {
1200// AliError(Form("Did not find eventType %s",eventType));
1201// delete bins;
1202// delete events;
1203// return 0x0;
1204// }
1205// delete events;
1206//
1207// Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
1208//
1209// if (ntrigger<=0)
1210// {
1211// AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
1212// delete bins;
1213// return 0x0;
1214// }
1215//
1216// AliAnalysisMuMuSpectra* sMinv = static_cast<AliAnalysisMuMuSpectra*>(spectraMinv.Clone());
1217// if(!sMinv)
1218// {
1219// AliError("Did not find inv mass spectra");
1220// }
1221//
1222//
1223// // binning.Print();
1224//
1225// AliAnalysisMuMuSpectra* spectra(0x0);
1226//
1227// AliAnalysisMuMuBinning::Range* bin;
1228// // AliAnalysisMuMuBinning::Range* binMinv;
1229// TIter next(bins);
1230//
1231// // TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
1232// // TIter nextFitType(fitTypeArray);
1233// // TObjString* fitType;
1234// TString flavour;
1235//
1236//
1237// while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1238// {
1239// std::cout << "---------------------------------" << std::endl;
1240// std::cout << "Fitting mean Pt in " << bin->AsString().Data() << " " << "for" << " " << eventType << "/" << trigger << "/" << centrality << "/" << pairCut << std::endl;
1241// std::cout << "---------------------------------" << std::endl;
1242//
1243//
1244// TString pname;
1245// if ( corrected ) pname = Form("MeanPtVsMinvUS%s_Corr",bin->AsString().Data());
1246// else pname = Form("MeanPtVsMinvUS%s",bin->AsString().Data());
1247//
1248// TProfile* hmPt = static_cast<TProfile*>(fMergeableCollection->GetObject(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),pname.Data()));
1249//
1250// if (!hmPt)
1251// {
1252// // if (!fBinning && bin->IsNullObject() )
1253// // {
1254// // // old file, we only had MinvUSPt
1255// // hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
1256// // }
1257//
1258// // if (!hmPt)
1259// // {
1260// AliDebug(1,Form("Could not find histo profile %s",pname.Data()));
1261// continue;
1262// // }
1263// }
1264// hmPt->Approximate();
1265//
1266//
1267// // hmPt = static_cast<TH1*>(hmPt->Clone(Form("mPtv%d",n++))); //Reuse this
1268//
1269// // if ( sbins.Contains("MULT") )
1270// // {
1271// // binMinv = binInt;
1272// // }
1273// // else
1274// // {
1275// // binMinv = bin;
1276// // }
1277//
1278// AliAnalysisMuMuJpsiResult* rMinv = static_cast<AliAnalysisMuMuJpsiResult*>(sMinv->GetResultForBin(*bin/*Minv*/));
1279//
1280// if( !rMinv )
1281// {
1282// AliError(Form("Did not find inv mass result inside spectra for bin %s",bin/*Minv*/->Quantity().Data()));
1283// }
1284//
1285// TObjArray* fits = rMinv->SubResults();
1286// AliAnalysisMuMuResult* fitMinv;
1287// TIter nextFit(fits);
1288// TString fitName;
1289//
1290// AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(*hmPt,
1291// trigger,
1292// eventType,
1293// pairCut,
1294// centrality,
1295// *bin);
1296//
1297// r->SetNofTriggers(ntrigger);
1298//
1299// nextFit.Reset();
1300//
1301// // new TCanvas;
1302// Int_t i = 1;
1303//
1304// while ( ( fitMinv = static_cast<AliAnalysisMuMuJpsiResult*>(nextFit())) )
1305// {
1306// std::cout << "" << std::endl;
1307// std::cout << "---------------" << "Fit " << i << "------------------" << std::endl;
1308// i++;
1309//
1310// fitName = fitMinv->Alias();
1311// std::cout << "" << std::endl;
1312// std::cout << "Getting signal parameters from " << eventType << "/" << trigger << "/" << centrality << "/" << pairCut << "/" << spectraMinv.GetName() << std::endl;
1313// std::cout << "" << std::endl;
1314//
1315// if(fitName.BeginsWith("JPSI2CB2VWG") || fitName.BeginsWith("MCTAILS"))
1316// {
1317// Double_t par[14] = {fitMinv->GetValue("kVWG"),fitMinv->GetValue("mVWG"),fitMinv->GetValue("sVWG1"),
1318// fitMinv->GetValue("sVWG2"),fitMinv->GetValue("kPsi"),fitMinv->GetValue("MeanJpsi"),fitMinv->GetValue("SigmaJpsi"),
1319// fitMinv->GetValue("alPsi"),fitMinv->GetValue("nlPsi"),fitMinv->GetValue("auPsi"),fitMinv->GetValue("nuPsi"),fitMinv->GetValue("kPsi'"),fitMinv->GetValue("NofJpsi"),fitMinv->GetErrorStat("NofJpsi")}; //Create a method to get the parameters //The last 2 parameters are not used in the fit
1320//
1321// if(fitName.BeginsWith("JPSI2CB2VWG"))
1322// {
1323// std::cout << "PRE-SET TAILS" << std::endl;
1324// std::cout << "" << std::endl;
1325// r->AddMeanPtFit("JPSI2CB2VWG:2",&par[0]); // FIXME: :2 should be in the default fit types but for meanpt fit naming is not there, change it
1326// }
1327// else
1328// {
1329// std::cout << "MC TAILS" << std::endl;
1330// std::cout << "" << std::endl;
1331// r->AddMeanPtFit("MCTAILS-JPSI2CB2VWG:2",&par[0]);
1332// }
1333//
1334//
1335// }
1336//
1337// //Make a part for NA48 and the rest of fitting functions.
1338//
1339// }
1340//
1341// flavour = bin->Flavour();
1342//
1343// if (!spectra)
1344// {
1345// TString spectraName(Form("MeanPtVsMinvUS-%s",binning.GetName()));
1346// if ( flavour.Length() > 0 )
1347// {
1348// spectraName += "-";
1349// spectraName += flavour;
1350// }
1351// if ( corrected )
1352// {
1353// spectraName += "-";
1354// spectraName += "Corr";
1355// }
1356//
1357// spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
1358// }
1359//
1360// spectra->AdoptResult(*bin,r);
1361//
1362// // if ( IsSimulation() )
1363// // {
1364// // SetNofInputParticles(*r);
1365// // }
1366//
1367//
1368// } // loop on bins
1369//
1370// // delete fitTypeArray;
1371// delete sMinv;
1372// delete bins;
1373//
1374// return spectra;
1375//}
1376
1377//_____________________________________________________________________________
1378void
1379AliAnalysisMuMu::GetParametersFromMC(TString& fitType, const char* pathCentrPairCut, const char* spectraName, AliAnalysisMuMuBinning::Range* bin) const
1380{
1381 /// Get the parameters from the associated simulation. Is intended to be used in minv fits, where we need the tails from JPsi (and Psi')
1382 // We can choose between the 2 associated simulations (a JPsi and Psi' ones) by setting the sim variable to "sim1" (fAssociatedSimulation by default) or "sim2" (fAssociatedSimulation2)
1383
1384 if ( !SIM() && !SIM2() )
1385 {
1386 AliError("Cannot get MC tails without associated simulation(s) file(s) !");
1387 fitType = "";
1388 return;
1389 }
1390
1391 TString subResultName("");
1392 if ( fitType.Contains("NA60NEW",TString::kIgnoreCase) ) subResultName = "PSINA60NEW";//_1.5_4.2
1393 else if ( fitType.Contains("CB2",TString::kIgnoreCase) ) subResultName = "PSICB2";//_2.2_3.9
1394 else
1395 {
1396 AliError("I don't know from which MC subresult take the tails");
1397 return;
1398 }
1399
1400 TObjArray* simArr = new TObjArray;
1401 if ( SIM() ) simArr->Add(SIM());
1402 if ( SIM2() && fitType.Contains("INDEPTAILS",TString::kIgnoreCase) ) simArr->Add(SIM2()); // If we have set the key to get the JPsi ans PsiP tails
1403
1404 TIter nextSim(simArr);
1405 AliAnalysisMuMu* currentSIM;
1406
1407 TString spath(pathCentrPairCut);
1408
1409 spath.Prepend(Form("/%s",First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kTRUE)).Data()));//FIXME: Care with this when there is more than one selection in the list
1410 spath.Prepend(Form("/%s",First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kTRUE)).Data()));
1411
1412
1413 while ( (currentSIM = static_cast<AliAnalysisMuMu*>(nextSim())) )
1414 {
1415 AliAnalysisMuMuSpectra* minvMCSpectra = currentSIM->SPECTRA(Form("%s/%s",spath.Data(),spectraName));
1416 if (!minvMCSpectra)
1417 {
1418 AliError(Form("Could not find spectra %s/%s for associated simulation",spath.Data(),spectraName));
1419 currentSIM->OC()->Print("*:Ali*");
1420 fitType = "";
1421 continue;
1422 }
1423
1424 AliAnalysisMuMuJpsiResult* minvMCResult = static_cast<AliAnalysisMuMuJpsiResult*>(minvMCSpectra->GetResultForBin(*bin));
1425 if ( !minvMCResult )
1426 {
1427 AliError(Form("Cannot get MC tails cause the minv result for bin %s in %s/%s was not found",bin->AsString().Data(),spath.Data(),spectraName));
1428 fitType = "";
1429 continue;
1430 }
1431
1432 AliAnalysisMuMuJpsiResult* r = dynamic_cast<AliAnalysisMuMuJpsiResult*>(minvMCResult->SubResult(subResultName.Data())); //FIXME: Find an independet way of naming results
1433 if ( r && subResultName.Contains("PSICB2") )
1434 {
1435 fitType += Form(":al%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("al%s",currentSIM->GetParticleName())));
1436 fitType += Form(":nl%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("nl%s",currentSIM->GetParticleName())));
1437 fitType += Form(":au%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("au%s",currentSIM->GetParticleName())));
1438 fitType += Form(":nu%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("nu%s",currentSIM->GetParticleName())));
1439
1440 std::cout << " Using MC " << currentSIM->GetParticleName() << " CB2 tails... " << std::endl;
1441 std::cout << std::endl;
1442 }
1443 else if ( r && subResultName.Contains("PSINA60NEW") )
1444 {
1445 fitType += Form(":p1L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p1L%s",currentSIM->GetParticleName())));
1446 fitType += Form(":p2L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p2L%s",currentSIM->GetParticleName())));
1447 fitType += Form(":p3L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p3L%s",currentSIM->GetParticleName())));
1448 fitType += Form(":p1R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p1R%s",currentSIM->GetParticleName())));
1449 fitType += Form(":p2R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p2R%s",currentSIM->GetParticleName())));
1450 fitType += Form(":p3R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p3R%s",currentSIM->GetParticleName())));
1451
1452 fitType += Form(":aL%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("aL%s",currentSIM->GetParticleName())));
1453 fitType += Form(":aR%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("aR%s",currentSIM->GetParticleName())));
1454
1455 std::cout << " Using MC " << currentSIM->GetParticleName() << " NA60New tails... " << std::endl;
1456 std::cout << std::endl;
1457 }
1458 else
1459 {
1460 AliError(Form("Cannot get MC tails. MC Subresult %s not found",minvMCResult->GetName()));
1461 fitType = "";
1462 }
1463 }
1464 delete simArr;
1465}
1466
1467//_____________________________________________________________________________
1468void
1469AliAnalysisMuMu::GetParametersFromResult(TString& fitType, AliAnalysisMuMuJpsiResult* minvResult) const
1470{
1471 // Gets the parameters from a result, is intended to be used for mean pt fits where we need the signal and backgroud parameters
1472
1473 AliWarning("Re-implement me !!!"); //FIXME: The parameters to get will depend on the fit function and also in this way is not suitable for other particles (ie Upsilon)(Find a way to get the particle(s) name)
1474
1475 TString msg("");
1476 if ( minvResult )
1477 {
1478 // We take the usual parameters (the ones from JPsi and the normalization of the Psi')
1479 fitType += Form(":kJPsi=%f",minvResult->GetValue("kJPsi")); //FIXME: Names are not correct
1480 fitType += Form(":mJPsi=%f",minvResult->GetValue("mJPsi"));
1481 fitType += Form(":sJPsi=%f",minvResult->GetValue("sJPsi"));
1482
1483 fitType += Form(":NofJPsi=%f",minvResult->GetValue("NofJPsi"));
1484 fitType += Form(":ErrStatNofJPsi=%f",minvResult->GetErrorStat("NofJPsi"));
1485
1486 fitType += Form(":kPsiP=%f",minvResult->GetValue("kPsiP"));
1487
1488 TString minvName(minvResult->GetName());
1489
1490 TString minvRangeParam = minvName;
1491 minvRangeParam.Remove(0,minvRangeParam.First("_") + 1);
1492 fitType += Form(":MinvRS=%s",minvRangeParam.Data());
1493
1494 fitType += Form(":FSigmaPsiP=%f",minvResult->GetValue("FSigmaPsiP"));
1495
1496 if ( fitType.Contains("CB2",TString::kIgnoreCase) )
1497 {
1498 fitType += Form(":alJPsi=%f",minvResult->GetValue("alJPsi"));
1499 fitType += Form(":nlJPsi=%f",minvResult->GetValue("nlJPsi"));
1500 fitType += Form(":auJPsi=%f",minvResult->GetValue("auJPsi"));
1501 fitType += Form(":nuJPsi=%f",minvResult->GetValue("nuJPsi"));
1502
1503 msg += "JPsi CB2 signal parameters";
1504 // fitType += Form(":NofPsiP=%f",minvResult->GetValue("NofPsiP"));
1505 // fitType += Form(":ErrStatNofPsiP=%f",minvResult->GetErrorStat("NofPsiP"));
1506
1507 if ( fitType.Contains("INDEPTAILS") )
1508 {
1509// minvName = minvResult->GetName();
1510 if ( minvName.Contains("INDEPTAILS") )
1511 {
1512 // In case we use independent parameters tails for JPsi and Psi' we take also the Psi' ones
1513 fitType += Form(":alPsiP=%f",minvResult->GetValue("alPsiP"));
1514 fitType += Form(":nlPsiP=%f",minvResult->GetValue("nlPsiP"));
1515 fitType += Form(":auPsiP=%f",minvResult->GetValue("auPsiP"));
1516 fitType += Form(":nuPsiP=%f",minvResult->GetValue("nuPsiP"));
1517 fitType += Form(":mPsiP=%f",minvResult->GetValue("mPsiP"));
1518 fitType += Form(":sPsiP=%f",minvResult->GetValue("sPsiP"));
1519
1520 msg += " + PsiP CB2 signal parameters";
1521 }
1522 else
1523 {
1524 AliError(Form("Cannot get PsiP tails from result. Result %s does not contain PsiP tails info => Fit will fail ",minvResult->GetName()));
1525 fitType = "";
1526 return;
1527 }
1528 }
1529 }
1530 else if ( fitType.Contains("NA60NEW",TString::kIgnoreCase) )
1531 {
1532 fitType += Form(":p1LJPsi=%f",minvResult->GetValue("p1LJPsi"));
1533 fitType += Form(":p2LJPsi=%f",minvResult->GetValue("p2LJPsi"));
1534 fitType += Form(":p3LJPsi=%f",minvResult->GetValue("p3LJPsi"));
1535 fitType += Form(":p1RJPsi=%f",minvResult->GetValue("p1RJPsi"));
1536 fitType += Form(":p2RJPsi=%f",minvResult->GetValue("p2RJPsi"));
1537 fitType += Form(":p3RJPsi=%f",minvResult->GetValue("p3RJPsi"));
1538
1539 fitType += Form(":aLJPsi=%f",minvResult->GetValue("aLJPsi"));
1540 fitType += Form(":aRJPsi=%f",minvResult->GetValue("aRJPsi"));
1541
1542 msg += "JPsi NA60New signal parameters";
1543
1544 if ( fitType.Contains("INDEPTAILS") )
1545 {
1546// TString minvName(minvResult->GetName());
1547 if ( minvName.Contains("INDEPTAILS") )
1548 {
1549 // In case we use independent parameters tails for JPsi and Psi' we take also the Psi' ones
1550 fitType += Form(":p1LPsiP=%f",minvResult->GetValue("p1LPsiP"));
1551 fitType += Form(":p2LPsiP=%f",minvResult->GetValue("p2LPsiP"));
1552 fitType += Form(":p3LPsiP=%f",minvResult->GetValue("p3LPsiP"));
1553 fitType += Form(":p1RPsiP=%f",minvResult->GetValue("p1RPsiP"));
1554 fitType += Form(":p2RPsiP=%f",minvResult->GetValue("p2RPsiP"));
1555 fitType += Form(":p3RPsiP=%f",minvResult->GetValue("p3RPsiP"));
1556
1557 fitType += Form(":aLPsiP=%f",minvResult->GetValue("aLPsiP"));
1558 fitType += Form(":aRPsiP=%f",minvResult->GetValue("aRPsiP"));
1559
1560 msg += " + PsiP NA60New signal parameters";
1561
1562 }
1563 else
1564 {
1565 AliError(Form("Cannot get PsiP tails from result. Result %s does not contain PsiP tails info => Fit will fail ",minvResult->GetName()));
1566 fitType = "";
1567 return;
1568 }
1569 }
1570 }
1571 else
1572 {
1573 AliError(Form("Cannot get the parameters from %s",minvResult->GetName()));
1574 fitType = "";
1575 return;
1576 }
1577 // Now we take the background parameters
1578 if ( fitType.Contains("VWG_") || fitType.Contains("VWGINDEPTAILS") ) //FIXME: Check that cannot be misunderstood(like Exp x Pol2..). In fact it can be misunderstood since the meanpt function name has also the name of the function to fit the bkg (free parameters). Also add the rest of the BKG functions
1579 {
1580 fitType += Form(":kVWG=%f",minvResult->GetValue("kVWG"));
1581 fitType += Form(":mVWG=%f",minvResult->GetValue("mVWG"));
1582 fitType += Form(":sVWG1=%f",minvResult->GetValue("sVWG1"));
1583 fitType += Form(":sVWG2=%f",minvResult->GetValue("sVWG2"));
1584
1585 msg += " + VWG Bkg parameters";
1586 }
1587 else if ( fitType.Contains("POL2EXP_") || fitType.Contains("POL2EXPINDEPTAILS") )
1588 {
1589 fitType += Form(":kPol2Exp=%f",minvResult->GetValue("kPol2Exp"));
1590 fitType += Form(":pol0=%f",minvResult->GetValue("pol0"));
1591 fitType += Form(":pol1=%f",minvResult->GetValue("pol1"));
1592 fitType += Form(":pol2=%f",minvResult->GetValue("pol2"));
1593 fitType += Form(":exp=%f",minvResult->GetValue("exp"));
1594
1595 msg += " + Pol2xExp Bkg parameters";
1596 }
1597 else if ( fitType.Contains("POL4EXP_") || fitType.Contains("POL4EXPINDEPTAILS") )
1598 {
1599 fitType += Form(":kPol4Exp=%f",minvResult->GetValue("kPol4Exp"));
1600 fitType += Form(":pol0=%f",minvResult->GetValue("pol0"));
1601 fitType += Form(":pol1=%f",minvResult->GetValue("pol1"));
1602 fitType += Form(":pol2=%f",minvResult->GetValue("pol2"));
1603 fitType += Form(":pol3=%f",minvResult->GetValue("pol3"));
1604 fitType += Form(":pol4=%f",minvResult->GetValue("pol4"));
1605 fitType += Form(":exp=%f",minvResult->GetValue("exp"));
1606
1607 msg += " + Pol4xExp Bkg parameters";
1608 }
1609 std::cout << "Using " << msg.Data() << " from " << minvResult->GetName() << " inv mass result" << std::endl;
1610 std::cout << "" << std::endl;
1611 }
1612 else
1613 {
1614 AliError(Form("Cannot get tails from result. Result %s not found",minvResult->GetName()));
1615 fitType = "";
1616 return;
1617 }
1618}
1619
1620//_____________________________________________________________________________
1621ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
1622{
1623 // Get the expected (from OCDB scalers) trigger count
1624
1625 AliAnalysisTriggerScalers ts(runNumber,Config()->OCDBPath());
1626
1627 TObjArray* triggers = TString(triggerList).Tokenize(",");
1628 TObjString* trigger;
1629 TIter next(triggers);
1630 ULong64_t n(0);
1631
1632 while ( ( trigger = static_cast<TObjString*>(next()) ) )
1633 {
1634 AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
1635 if (item)
1636 {
1637 n += item->Value();
1638 }
1639 delete item;
1640 }
1641 delete triggers;
1642
1643 return n;
1644}
1645
1646//_____________________________________________________________________________
1647AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
1648{
1649 /// get a given spectra
1650
1651 TString swhat(what);
1652 TString sflavour(flavour);
1653 swhat.ToUpper();
1654 sflavour.ToUpper();
1655
1656 TString spectraName(Form("/%s/%s/PP/%s/PSI-%s",
1657 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation())).Data(),
1658 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation())).Data(),
1659 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation())).Data(),
1660 swhat.Data()));
1661
1662 if (sflavour.Length()>0)
1663 {
1664 spectraName += "-";
1665 spectraName += sflavour.Data();
1666 }
1667
1668 return SPECTRA(spectraName.Data());
1669}
1670
1671//_____________________________________________________________________________
1672TH1* AliAnalysisMuMu::PlotAccEfficiency(const char* whatever)
1673{
1674 //FIXME::Make it general
1675 if ( !IsSimulation() )
1676 {
1677 AliError("Could not get AccxEff histo: Not simulation file");
1678 return 0x0;
1679 }
1680
1681 TString path(Form("/%s/%s/%s/%s",
1682 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kTRUE)).Data(),
1683 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kTRUE)).Data(),
1684 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kTRUE)).Data(),
1685 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kTRUE)).Data()));
1686
1687 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("%s/%s",path.Data(),whatever)));
1688 if ( !s )
1689 {
1690 AliError(Form("No AccEff spectra %s found in %s",whatever,path.Data()));
1691 return 0x0;
1692 }
1693
1694 return s->Plot(Form("AccEff%s",GetParticleName()),"PSICOUNT",kFALSE);//_2.2_3.9
1695
1696}
1697
1698//_____________________________________________________________________________
1699TH1* AliAnalysisMuMu::PlotSystematicsTestsRelative(const char* quantity,const char* flavour,const char* value2Test)
1700{
1701 /// what,quantity and flavour defines de binning to test (output will be 1 plot per bin)
1702 /// value2test is the observable we want to test ( i.e. NJpsi(bin)/integratedNJpsi, <pt>(bin)/integrated<pt>... )
1703 /// --value2test == yield -> NJpsi(bin)/integratedNJpsi
1704 /// --value2test == mpt -> <pt>(bin)/integrated<pt>
1705
1706 TString path(Form("%s/%s/%s/%s",
1707 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
1708 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
1709 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
1710 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
1711
1712 TString svalue2Test(value2Test);
1713 TString shName("");
1714 TString sObsInt("");
1715 TString sObsBin("");
1716 TString sflavour(flavour);
1717 TString squantity(quantity);
1718 squantity.ToUpper();
1719
1720 if ( !svalue2Test.CompareTo("yield",TString::kIgnoreCase) )
1721 {
1722 sObsInt = "PSI-INTEGRATED-AccEffCorr";
1723 sObsBin = Form("PSI-%s-AccEffCorr",squantity.Data());
1724 svalue2Test = "NofJPsi";
1725 shName = "N^{J/#psi}_{bin}/N^{J/#psi}_{int}";
1726 }
1727 else if ( !svalue2Test.CompareTo("mpt",TString::kIgnoreCase) )
1728 {
1729 sObsInt = "PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS";
1730 sObsBin = Form("PSI-%s-AccEffCorr-MeanPtVsMinvUS",squantity.Data());
1731 svalue2Test = "MeanPtJPsi";
1732 shName = "<p_{t}>^{bin}/<p_{t}>^{int}";
1733 }
1734 else
1735 {
1736 AliError("unrecognized value to test");
1737 return 0x0;
1738 }
1739
1740 TString id(Form("/TESTSYST/%s",path.Data()));
1741
1742 //--Get the integrated results
1743 AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),sObsInt.Data())));
1744 if ( !sInt )
1745 {
1746 AliError(Form("No integrated spectra %s found in %s",sObsInt.Data(),path.Data()));
1747 return 0x0;
1748 }
1749
1750 TObjArray* bin = BIN()->CreateBinObjArray("psi","integrated","");
1751 if ( !bin )
1752 {
1753 AliError("No integrated bin found");
1754 return 0x0;
1755 }
1756 AliAnalysisMuMuBinning::Range* r = static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(0));
1757
1758 AliAnalysisMuMuJpsiResult* resInt = static_cast<AliAnalysisMuMuJpsiResult*>(sInt->GetResultForBin(*r));
1759 if ( !resInt )
1760 {
1761 AliError(Form("No integrated result found in spectra %s at %s",sObsInt.Data(),path.Data()));
1762 return 0x0;
1763 }
1764 TObjArray* sresIntArray = resInt->SubResults();
1765// Int_t nsresInt = sresIntArray->GetEntriesFast();
1766
1767 bin = BIN()->CreateBinObjArray("psi",squantity.Data(),sflavour.Data());//memory leak?
1768 if ( !bin )
1769 {
1770 AliError(Form("%s-%s-%s binning does not exist","psi",squantity.Data(),sflavour.Data()));
1771 return 0x0;
1772 }
1773
1774 Int_t nbin = bin->GetEntries();
1775 Int_t nsres = sresIntArray->GetEntries();
1776 TObjArray* sResultNameArray= new TObjArray();
1777 std::vector<std::vector<double> > valuesArr;
1778 valuesArr.resize(nbin+1, std::vector<double>(nsres,0));
1779 std::vector<std::vector<double> > valuesErrorArr;
1780 valuesErrorArr.resize(nbin+1, std::vector<double>(nsres,0));
1781
1782 TIter nextIntSubResult(sresIntArray);
1783 AliAnalysisMuMuResult* sresInt(0x0);
1784 Int_t nBkgParametrizations(0);
1785 Int_t j(0);
1786 while ( ( sresInt = static_cast<AliAnalysisMuMuResult*>(nextIntSubResult()) ) ) //Integrated SubResult loop
1787 {
1788 TString sresIntName(sresInt->GetName());
1789
1790 valuesArr[0][j] = sresInt->GetValue(svalue2Test.Data());
1791 valuesErrorArr[0][j] = sresInt->GetErrorStat(svalue2Test.Data());
1792 sResultNameArray->Add(new TObjString(sresIntName.Data()));
1793
1794 if ( sresIntName.Contains("CB2") )
1795 {
1796 nBkgParametrizations++;
1797 }
1798 j++;
1799 }
1800
1801 //--Get the bin per bin results
1802 AliAnalysisMuMuSpectra* sBin = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),sObsBin.Data())));
1803 if ( !sBin )
1804 {
1805 AliError(Form("No integrated spectra %s found in %s",sObsBin.Data(),path.Data()));
1806 return 0x0;
1807 delete bin;
1808 }
1809
1810 TIter nextBin(bin);
1811 AliAnalysisMuMuJpsiResult* sresBin(0x0);
1812 Int_t i(1);
1813 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
1814 {
1815 sresBin = static_cast<AliAnalysisMuMuJpsiResult*>(sBin->GetResultForBin(*r));
1816 if ( !sresBin )
1817 {
1818 AliError(Form("No result found in spectra %s at %s for bin %s",sObsInt.Data(),path.Data(),r->AsString().Data()));
1819 return 0x0;
1820 delete bin;
1821 }
1822
1823 TObjArray* sresBinArray = sresBin->SubResults();
1824
1825 TIter nextSubResult(sresBinArray);
1826 j = 0;
1827 while ( (sresBin = static_cast<AliAnalysisMuMuJpsiResult*>(nextSubResult())) ) // Subresults loop
1828 {
1829 valuesArr[i][j] = sresBin->GetValue(svalue2Test.Data());
1830 valuesErrorArr[i][j] = sresBin->GetErrorStat(svalue2Test.Data());
1831
1832 j++;
1833
1834 } //End Subresults loop
1835 i++;
1836
1837 } //End bin loop
1838
1839 TH1* hsyst = new TH1F(Form("%s_Systematics",value2Test),Form("%s Systematics results",value2Test),nbin,0,nbin);
1840
1841 TString binName("");
1842 for ( Int_t b = 1 ; b < i ; b++ ) //Bin loop
1843 {
1844 binName = static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(b-1))->AsString().Data();
1845 TString savePath(Form("%s/%s",id.Data(),binName.Data()));
1846
1847 Int_t binUCTotal(1);
1848
1849 TH1* hratiosBin = new TH1D(Form("SystTests_%s_Bkg_%s",sObsBin.Data(),binName.Data()),
1850 Form("%s Systematics tests for %s",binName.Data(),shName.Data()),j*nBkgParametrizations,0,
1851 j*nBkgParametrizations);
1852
1853 for ( Int_t k = 0 ; k <= j ; k++ )
1854 {
1855 TString hName("");
1856
1857 if ( k == 0 ) hName = "UC";
1858 else hName = sResultNameArray->At(k-1)->GetName();
1859
1860
1861 TH1* hratios = new TH1D(Form("SystTests_%s_%s",sObsBin.Data(),hName.Data()),
1862 Form("%s Systematics tests for %s(%s)",binName.Data(),shName.Data(),hName.Data()),j,0,j);
1863 hratios->SetNdivisions(j+1);
1864
1865 TH1* hratiosUC(0x0);
1866 if ( k != 0 )
1867 {
1868 hratiosUC = new TH1D(Form("SystTests_%s_%s_Bkg",sObsBin.Data(),hName.Data()),
1869 Form("%s Systematics tests for %s(%s bkg)",binName.Data(),shName.Data(),hName.Data()),nBkgParametrizations,0,nBkgParametrizations);
1870 hratiosUC->SetNdivisions(nBkgParametrizations+1);
1871 }
1872
1873 TString signalName(hName.Data());
1874 Int_t sizeName = signalName.Sizeof();
1875 signalName.Remove(2,sizeName-3);
1876 Int_t binUC(1);
1877
1878 for ( Int_t l = 0; l < j ; l++) //Subresults loop
1879 {
1880 TString binSignalName(sResultNameArray->At(l)->GetName());
1881
1882 Double_t ratio,ratioError;
1883 if ( k == 0)
1884 {
1885 ratio = valuesArr[b][l] / valuesArr[0][l];
1886 ratioError = TMath::Sqrt( TMath::Power(valuesErrorArr[b][l] / valuesArr[0][l],2.) + TMath::Power(valuesArr[b][l]*valuesErrorArr[0][l] / TMath::Power(valuesArr[0][l],2.),2.) );
1887
1888 hratios->GetXaxis()->SetBinLabel(l+1,sResultNameArray->At(l)->GetName());
1889 }
1890 else
1891 {
1892 ratio = valuesArr[b][l] / valuesArr[0][k-1];
1893 ratioError = TMath::Sqrt( TMath::Power(valuesErrorArr[b][l] / valuesArr[0][k-1],2.) + TMath::Power(valuesArr[b][l]*valuesErrorArr[0][k-1] / TMath::Power(valuesArr[0][k-1],2.),2.) );
1894
1895 hratios->GetXaxis()->SetBinLabel(l+1,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1896
1897 if ( binSignalName.Contains(signalName.Data()) )
1898 {
1899 hratiosUC->GetXaxis()->SetBinLabel(binUC,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1900
1901 hratiosUC->SetBinContent(binUC,ratio);
1902 hratiosUC->SetBinError(binUC,ratioError);
1903
1904 hratiosBin->GetXaxis()->SetBinLabel(binUCTotal,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1905
1906 hratiosBin->SetBinContent(binUCTotal,ratio);
1907 hratiosBin->SetBinError(binUCTotal,ratioError);
1908
1909 binUC++;
1910 binUCTotal++;
1911 }
1912 }
1913
1914 hratios->SetBinContent(l+1,ratio);
1915 hratios->SetBinError(l+1,ratioError);
1916 }
1917
1918
1919
1920
1921// TH1* o = OC()->Histo(Form("%s",savePath.Data()),hratios->GetName());
1922//
1923// if (o)
1924// {
1925// AliWarning(Form("Replacing %s/%s",savePath.Data(),hratios->GetName()));
1926// OC()->Remove(Form("%s/%s",savePath.Data(),hratios->GetName()));
1927// }
1928//
1929// Bool_t adoptOK = OC()->Adopt(savePath.Data(),hratios);
1930//
1931// if ( adoptOK ) std::cout << "+++syst histo " << hratios->GetName() << " adopted" << std::endl;
1932// else AliError(Form("Could not adopt syst histo %s",hratios->GetName()));
1933//
1934// if ( hratiosUC )
1935// {
1936// o = OC()->Histo(Form("%s",savePath.Data()),hratiosUC->GetName());
1937//
1938// if (o)
1939// {
1940// AliWarning(Form("Replacing %s/%s",savePath.Data(),hratiosUC->GetName()));
1941// OC()->Remove(Form("%s/%s",savePath.Data(),hratiosUC->GetName()));
1942// }
1943//
1944// adoptOK = OC()->Adopt(savePath.Data(),hratiosUC);
1945//
1946// if ( adoptOK ) std::cout << "+++syst histo " << hratiosUC->GetName() << " adopted" << std::endl;
1947// else AliError(Form("Could not adopt syst histo %s",hratiosUC->GetName()));
1948// }
1949 }
1950
1951
1952 //Syst computation for bin
1953 Double_t num(0.),deno(0.);
1954 for ( Int_t m = 1 ; m <= hratiosBin->GetNbinsX() ; m++ )
1955 {
1956 Double_t value = hratiosBin->GetBinContent(m);
1957 Double_t error = hratiosBin->GetBinError(m);
1958
1959 num += value/TMath::Power(error,2.);
1960 deno += 1./TMath::Power(error,2.);
1961 }
1962
1963 Double_t mean = num/deno;
1964 Double_t v1(0.),v2(0.),sum(0.);
1965 for ( Int_t l = 1 ; l <= hratiosBin->GetNbinsX() ; l++ )
1966 {
1967 Double_t value = hratiosBin->GetBinContent(l);
1968 Double_t error = hratiosBin->GetBinError(l);
1969
1970 Double_t wi = 1./TMath::Power(error,2.);
1971 v1 += wi;
1972 v2 += wi*wi;
1973 Double_t diff = value - mean;
1974 sum += wi*diff*diff;
1975
1976 }
1977
1978 Double_t syst = TMath::Sqrt( (v1/(v1*v1-v2)) * sum);
1979
1980 hsyst->GetXaxis()->SetBinLabel(b,binName.Data());
1981 hsyst->SetBinContent(b,(syst*100.)/mean);
1982
1983
1984 //___
1985 TF1* meanF = new TF1("mean","[0]",0,j*nBkgParametrizations);
1986 meanF->SetParameter(0,mean);
1987
1988 TF1* meanFPS = new TF1("meanPS","[0]",0,j*nBkgParametrizations);
1989 meanFPS->SetParameter(0,mean+syst);
1990 meanFPS->SetLineStyle(2);
1991
1992 TF1* meanFMS = new TF1("meanMS","[0]",0,j*nBkgParametrizations);
1993 meanFMS->SetParameter(0,mean-syst);
1994 meanFMS->SetLineStyle(2);
1995
1996 hratiosBin->GetListOfFunctions()->Add(meanF);
1997 hratiosBin->GetListOfFunctions()->Add(meanFPS);
1998 hratiosBin->GetListOfFunctions()->Add(meanFMS);
1999
2000 TH1* o = OC()->Histo(Form("%s",savePath.Data()),hratiosBin->GetName());
2001
2002 if (o)
2003 {
2004 AliWarning(Form("Replacing %s/%s",savePath.Data(),hratiosBin->GetName()));
2005 OC()->Remove(Form("%s/%s",savePath.Data(),hratiosBin->GetName()));
2006 }
2007
2008 Bool_t adoptOK = OC()->Adopt(savePath.Data(),hratiosBin);
2009
2010 if ( adoptOK ) std::cout << "+++syst histo " << hratiosBin->GetName() << " adopted" << std::endl;
2011 else AliError(Form("Could not adopt syst histo %s",hratiosBin->GetName()));
2012 }
2013
2014
2015 TH1* o = OC()->Histo(Form("%s",id.Data()),hsyst->GetName());
2016
2017 if (o)
2018 {
2019 AliWarning(Form("Replacing %s/%s",id.Data(),hsyst->GetName()));
2020 OC()->Remove(Form("%s/%s",id.Data(),hsyst->GetName()));
2021 }
2022
2023 Bool_t adoptOK = OC()->Adopt(id.Data(),hsyst);
2024
2025 if ( adoptOK ) std::cout << "+++syst histo " << hsyst->GetName() << " adopted" << std::endl;
2026 else AliError(Form("Could not adopt syst histo %s",hsyst->GetName()));
2027
2028
2029 delete bin;
2030 delete sResultNameArray;
2031
2032 return 0x0;
2033
2034}
2035
2036
2037//_____________________________________________________________________________
2038TH1* AliAnalysisMuMu::PlotJpsiYield(const char* whatever)
2039{
2040
2041 //FIXME::Make it general
2042 if ( IsSimulation() )
2043 {
2044 AliError("Cannot compute J/Psi yield: Is a simulation file");
2045 return 0x0;
2046 }
2047
2048 TString path(Form("/%s/%s/%s/%s",
2049 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
2050 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
2051 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
2052 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
2053
2054 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("%s/%s",path.Data(),whatever)));
2055 if ( !s )
2056 {
2057 AliError(Form("No spectra %s found in %s",whatever,path.Data()));
2058 return 0x0;
2059 }
2060
2061 std::cout << "Number of J/Psi:" << std::endl;
2062 TH1* hry = s->Plot("NofJPsi","PSIPSIPRIMECB2VWGINDEPTAILS",kFALSE); //Number of Jpsi
2063 std::cout << "" << std::endl;
2064
2065 std::cout << "Equivalent number of MB events:" << std::endl;
2066 TH1* hN = ComputeEquNofMB();
2067 std::cout << "" << std::endl;
2068
2069 TH1* hy = static_cast<TH1*>(hry->Clone("CorrJPsiYields"));
2070 Double_t bR = 0.0593; // BR(JPsi->mu+mu-)
2071 Double_t bRerror = 0.0006 ;
2072
2073 for (Int_t i = 1 ; i <= hy->GetNbinsX() ; i++)
2074 {
2075 Double_t yield = hry->GetBinContent(i)/(hN->GetBinContent(i)*bR);
2076 Double_t yieldError = TMath::Sqrt(TMath::Power(hry->GetBinError(i)/(hN->GetBinContent(i)*bR),2.) +
2077 TMath::Power(hN->GetBinError(i)*bR/TMath::Power(hN->GetBinContent(i)*bR,2.),2.) +
2078 TMath::Power(hry->GetBinContent(i)*hN->GetBinContent(i)*bRerror/TMath::Power(hN->GetBinContent(i)*bR,2.),2.));
2079
2080 std::cout << yield << " +- " << yieldError << std::endl;
2081
2082 hy->SetBinContent(i,yield);
2083 hy->SetBinError(i,yieldError);
2084 }
2085
2086 delete hry;
2087 delete hN;
2088
2089 return hy;
2090}
2091
2092
2093//_____________________________________________________________________________
2094UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
2095 const char* eventSelection, Int_t runNumber)
2096{
2097 TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
2098 TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
2099 TObjArray* a = TString(triggerList).Tokenize(" ");
2100 TIter next(a);
2101 TObjString* str;
2102
2103 UInt_t n(0);
2104
2105 TString sEventSelection(eventSelection);
2106 sEventSelection.ToUpper();
2107
2108 if ( kevent->FindObject(sEventSelection.Data()) )
2109 {
2110 while ( ( str = static_cast<TObjString*>(next()) ) )
2111 {
2112 if ( ktrigger->FindObject(str->String().Data()) )
2113 {
2114 if ( runNumber < 0 )
2115 {
2116 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
2117 }
2118 else
2119 {
2120 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
2121 }
2122 }
2123 }
2124 }
2125
2126 delete a;
2127 delete ktrigger;
2128 delete kevent;
2129 return n;
2130}
2131
2132//_____________________________________________________________________________
2133Bool_t
2134AliAnalysisMuMu::GetCollections(const char* rootfile,
2135 AliMergeableCollection*& oc,
2136 AliCounterCollection*& cc,
2137 AliAnalysisMuMuBinning*& bin,
2138 std::set<int>& runnumbers)
2139{
2140 oc = 0x0;
2141 cc = 0x0;
2142 bin = 0x0;
2143
2144 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
2145 if (!f)
2146 {
2147 f = TFile::Open(rootfile);
2148 }
2149
2150 if ( !f || f->IsZombie() )
2151 {
2152 return kFALSE;
2153 }
2154
2155 f->GetObject("OC",oc);
2156 if (!oc)
2157 {
2158 f->GetObject("MC",oc);
2159 }
2160 f->GetObject("CC",cc);
2161
2162 TIter next(f->GetListOfKeys());
2163 TKey* key;
2164
2165 while ( ( key = static_cast<TKey*>(next())) && !bin )
2166 {
2167 if ( strcmp(key->GetClassName(),"AliAnalysisMuMuBinning")==0 )
2168 {
2169 bin = dynamic_cast<AliAnalysisMuMuBinning*>(key->ReadObj());
2170 }
2171 }
2172
2173 if (!oc || !cc)
2174 {
2175 AliErrorClass("Old file. Please upgrade it!");
2176
2177 return kFALSE;
2178 }
2179
2180 // get run list
2181 TObjArray* runs = cc->GetKeyWords("run").Tokenize(",");
2182 runs->Sort();
2183 TIter nextRun(runs);
2184 TObjString* srun;
2185
2186 runnumbers.clear();
2187
2188 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2189 {
2190 runnumbers.insert(srun->String().Atoi());
2191 }
2192
2193 delete runs;
2194
2195 return kTRUE;
2196}
2197
2198//_____________________________________________________________________________
2199Bool_t AliAnalysisMuMu::IsSimulation() const
2200{
2201 // whether or not we have MC information
2202
2203// return kFALSE;
2204
2205 if (!fMergeableCollection) return kFALSE;
2206
2207 TList* list = fMergeableCollection->CreateListOfKeys(0);
2208 TIter next(list);
2209 TObjString* str;
2210 Bool_t ok(kFALSE);
2211
2212 while ( ( str = static_cast<TObjString*>(next()) ) )
2213 {
2214 if ( str->String().Contains(AliAnalysisMuMuBase::MCInputPrefix()) ) ok = kTRUE;
2215 }
2216 delete list;
2217
2218 return ok;
2219}
1afce1ce 2220
cb690a12 2221//_____________________________________________________________________________
2222Int_t
2223AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour, Bool_t fitmPt)
2224{
2225 // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
2226 // what="integrated" => fit only fully integrated MinvUS
2227 // what="pt" => fit MinvUS in pt bins
2228 // what="y" => fit MinvUS in y bins
2229 // what="pt,y" => fit MinvUS in (pt,y) bins
2230
2231 TStopwatch timer;
2232
2233 if (!fMergeableCollection)
2234 {
2235 AliError("No mergeable collection. Consider Upgrade()");
2236 return 0;
2237 }
2238
2239 Int_t nfits(0);
2240
2241 TObjArray* triggerArray = Config()->GetListElements(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation());
2242 TObjArray* eventTypeArray = Config()->GetListElements(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation());
2243 TObjArray* pairCutArray = Config()->GetListElements(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation());
2244 TObjArray* whatArray = TString(what).Tokenize(",");
2245 TObjArray* centralityArray = Config()->GetListElements(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation());
2246
2247 TIter nextTrigger(triggerArray);
2248 TIter nextEventType(eventTypeArray);
2249 TIter nextPairCut(pairCutArray);
2250 TIter nextWhat(whatArray);
2251 TIter nextCentrality(centralityArray);
2252
2253 TObjString* trigger;
2254 TObjString* eventType;
2255 TObjString* pairCut;
2256 TObjString* swhat;
2257 TObjString* centrality;
2258
2259 while ( ( swhat = static_cast<TObjString*>(nextWhat()) ) )
2260 {
2261 AliAnalysisMuMuBinning* binning(0x0);
2262
2263 if ( fBinning && swhat->String().Length() > 0 )
2264 {
2265 binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
2266 }
2267 else
2268 {
2269 binning = new AliAnalysisMuMuBinning;
2270 binning->AddBin("psi",swhat->String().Data());
2271 }
2272
2273 StdoutToAliDebug(1,std::cout << "++++++++++++ swhat=" << swhat->String().Data() << std::endl;);
2274
2275 std::cout << "" << std::endl;
2276 std::cout << "++++++++++++++++++" << "NEW BIN TYPE" << "+++++++++++++++++++" << std::endl;
2277 std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2278 std::cout << "++++++++++++ swhat=" << swhat->String().Data() << "++++++++++++++++++++" << std::endl;
2279 std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2280 std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2281 std::cout << "" << std::endl;
2282 std::cout << "" << std::endl;
2283
2284 if (!binning)
2285 {
2286 AliError("oups. binning is NULL");
2287 continue;
2288 }
2289
2290 StdoutToAliDebug(1,binning->Print(););
2291
2292 nextTrigger.Reset();
2293
2294 while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
2295 {
2296 AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
2297
2298 nextEventType.Reset();
2299
2300 while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
2301 {
2302 AliDebug(1,Form("--EVENTTYPE %s",eventType->String().Data()));
2303
2304 nextPairCut.Reset();
2305
2306 while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
2307 {
2308 AliDebug(1,Form("----PAIRCUT %s",pairCut->String().Data()));
2309
2310 nextCentrality.Reset();
2311
2312 while ( ( centrality = static_cast<TObjString*>(nextCentrality()) ) )
2313 {
2314 AliDebug(1,"----Fitting...");
2315
2316 AliAnalysisMuMuSpectra* spectra = FitParticle("psi",
2317 trigger->String().Data(), //Uncomment
2318 eventType->String().Data(),
2319 pairCut->String().Data(),
2320 centrality->String().Data(),
2321 *binning);
2322
2323 AliDebug(1,Form("----fitting done spectra = %p",spectra));
2324
2325 TString id(Form("/%s/%s/%s/%s",eventType->String().Data(),
2326 trigger->String().Data(),
2327 centrality->String().Data(),
2328 pairCut->String().Data()));
2329 TObject* o;
2330
2331 if ( spectra )
2332 {
2333 ++nfits;
2334
2335 o = fMergeableCollection->GetObject(id.Data(),spectra->GetName());
2336
2337 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2338
2339 if (o)
2340 {
2341 AliWarning(Form("Replacing %s/%s",id.Data(),spectra->GetName()));
2342 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectra->GetName()));
2343 }
2344
2345 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectra);
2346
2347 if ( adoptOK ) std::cout << "+++Spectra " << spectra->GetName() << " adopted" << std::endl;
2348 else AliError(Form("Could not adopt spectra %s",spectra->GetName()));
2349
2350 StdoutToAliDebug(1,spectra->Print(););
2351 }
2352 else AliError("Error creating spectra"); //Uncomment
2353
2354 AliDebug(1,"----Fitting corrected spectra...");
2355
2356 AliAnalysisMuMuSpectra* spectraCorr = FitParticle("psi",
2357 trigger->String().Data(),
2358 eventType->String().Data(),
2359 pairCut->String().Data(),
2360 centrality->String().Data(),
2361 *binning,"minv",kTRUE);
2362
2363 AliDebug(1,Form("----fitting done corrected spectra = %p",spectraCorr));
2364
2365 o = 0x0;
2366 if ( spectraCorr )
2367 {
2368 ++nfits;
2369
2370 o = fMergeableCollection->GetObject(id.Data(),spectraCorr->GetName());
2371
2372 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2373
2374 if (o)
2375 {
2376 AliWarning(Form("Replacing %s/%s",id.Data(),spectraCorr->GetName()));
2377 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraCorr->GetName()));
2378 }
2379
2380 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraCorr);
2381
2382 if ( adoptOK ) std::cout << "+++Spectra " << spectraCorr->GetName() << " adopted" << std::endl;
2383 else AliError(Form("Could not adopt spectra %s",spectraCorr->GetName()));
2384
2385 StdoutToAliDebug(1,spectraCorr->Print(););
2386 }
2387 else AliError("Error creating spectra");
2388
2389
2390 if (fitmPt)
2391 {
2392 AliDebug(1,"----Fitting mean pt...");
2393
2394 std::cout << "" << std::endl;
2395 std::cout << "" << std::endl;
2396 std::cout << "++++++++++++ Fitting mean Pt for " << swhat->String().Data() << " " << "slices" << std::endl; //Uncomment
2397
2398 if ( spectra )
2399 {
2400 AliAnalysisMuMuSpectra* spectraMeanPt = FitParticle("psi",
2401 trigger->String().Data(),
2402 eventType->String().Data(),
2403 pairCut->String().Data(),
2404 centrality->String().Data(),
2405 *binning,"mpt"/*,*spectra*/);
2406
2407
2408
2409 AliDebug(1,Form("----fitting done spectra = %p",spectraMeanPt));
2410 o = 0x0;
2411
2412 if ( spectraMeanPt )
2413 {
2414 ++nfits; //Review this
2415
2416 o = fMergeableCollection->GetObject(id.Data(),spectraMeanPt->GetName());
2417
2418 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2419
2420 if (o)
2421 {
2422 AliWarning(Form("Replacing %s/%s",id.Data(),spectraMeanPt->GetName()));
2423 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraMeanPt->GetName()));
2424 }
2425
2426 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraMeanPt);
2427 // spectraMeanPt->Print();
2428 if ( adoptOK ) std::cout << "+++Spectra " << spectraMeanPt->GetName() << " adopted" << std::endl;
2429 else AliError(Form("Could not adopt spectra %s",spectraMeanPt->GetName()));
2430 }
2431 else AliError("Error creating spectra");
2432
2433 }
2434 else std::cout << "Mean pt fit failed: No inv mass spectra for " << swhat->String().Data() << " " << "slices" << std::endl; //Uncomment
2435
2436
2437 std::cout << "++++++++++++ Fitting corrected mean Pt for" << " " << swhat->String().Data() << " " << "slices" << std::endl;
2438
2439 if ( spectraCorr )
2440 {
2441 AliAnalysisMuMuSpectra* spectraMeanPtCorr = FitParticle("psi",
2442 trigger->String().Data(),
2443 eventType->String().Data(),
2444 pairCut->String().Data(),
2445 centrality->String().Data(),
2446 *binning,"mpt"/*,*spectraCorr*/,kTRUE);
2447
2448
2449
2450 AliDebug(1,Form("----fitting done spectra = %p",spectraMeanPtCorr));
2451
2452 o = 0x0;
2453
2454 if ( spectraMeanPtCorr )
2455 {
2456 ++nfits; //Review this
2457
2458 o = fMergeableCollection->GetObject(id.Data(),spectraMeanPtCorr->GetName());
2459
2460 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2461
2462 if (o)
2463 {
2464 AliWarning(Form("Replacing %s/%s",id.Data(),spectraMeanPtCorr->GetName()));
2465 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraMeanPtCorr->GetName()));
2466 }
2467
2468 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraMeanPtCorr);
2469 // spectraMeanPtCorr->Print();
2470 if ( adoptOK ) std::cout << "+++Spectra " << spectraMeanPtCorr->GetName() << " adopted" << std::endl;
2471 else AliError(Form("Could not adopt spectra %s",spectraMeanPtCorr->GetName()));
2472
2473 }
2474 else AliError("Error creating spectra");
2475
2476 }
2477
2478 else std::cout << "Corrected mean pt fit failed: No corrected inv mass spectra for " << swhat->String().Data() << " " << "slices" << std::endl;
2479
2480 }
2481 }
2482 }
2483 }
2484 }
2485 }
2486
2487 delete whatArray;
2488 delete triggerArray;
2489 delete eventTypeArray;
2490 delete pairCutArray;
2491 delete centralityArray;
2492
2493 StdoutToAliDebug(1,timer.Print(););
1afce1ce 2494
cb690a12 2495 if (nfits)
2496 {
2497 Update();
2498// ReOpen(fFilename,"UPDATE");
2499// fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
2500// ReOpen(fFilename,"READ");
2501 }
2502
2503
2504 return nfits;
2505
2506}
2507
2508//_____________________________________________________________________________
2509TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const char* event1,
2510 const char* trigger2, const char* event2,
2511 Bool_t drawFills,
2512 Bool_t asRejection) const
2513{
2514 if (!CC()) return 0x0;
2515
2516 const std::set<int>& runnumbers = RunNumbers();
2517
2518 TGraphErrors* g = new TGraphErrors(runnumbers.size());
2519
2520 std::set<int>::const_iterator it;
2521 Int_t i(0);
2522
2523 Double_t ymin(TMath::Limits<double>::Max());
2524 Double_t ymax(TMath::Limits<double>::Min());
2525
2526 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2527 {
2528 Int_t runNumber = *it;
2529 Double_t n = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger1,event1,runNumber));
2530 Double_t d = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger2,event2,runNumber));
2531 if (n>0 && d>0)
2532 {
2533 Double_t y = n/d;
2534
2535 if ( fCorrectionPerRun )
2536 {
2537 Double_t xcorr,ycorr;
2538 fCorrectionPerRun->GetPoint(i,xcorr,ycorr); // note that the fact that xcorr==runNumber has been checked by the SetCorrectionPerRun method
2539 y *= ycorr;
2540 // FIXME: should get the correction error here
2541 }
2542
2543 if ( asRejection ) y = 100*(1.0 - y);
2544 ymin = TMath::Min(ymin,y);
2545 ymax = TMath::Max(ymax,y);
2546 Double_t yerr = y*AliAnalysisMuMuResult::ErrorAB(n,TMath::Sqrt(n),d,TMath::Sqrt(d));
2547 g->SetPoint(i,runNumber,y);
2548 g->SetPointError(i,0.5,yerr);
2549
2550 ++i;
2551 }
2552
2553 }
2554
2555 TH2* hframe = new TH2F(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2556 Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2557 runnumbers.size()+50,
2558 *(runnumbers.begin())-25,
2559 *(runnumbers.rbegin())+25,100,0,ymax*1.3);
2560
2561 gStyle->SetOptStat(0);
2562
2563 hframe->Draw();
2564
2565 hframe->GetXaxis()->SetNoExponent();
2566
2567 hframe->GetYaxis()->SetTitle(asRejection ? "Rejection (%)" : "Ratio");
2568
2569 g->Set(i);
2570 g->SetTitle(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2));
2571 g->GetXaxis()->SetNoExponent();
2572 g->Draw("lp");
2573
2574 AliAnalysisTriggerScalers ts(RunNumbers(),Config()->OCDBPath());
2575
2576 if ( drawFills )
2577 {
2578 ts.DrawFills(ymin,ymax);
2579 g->Draw("lp");
2580 }
2581
2582
2583 std::map<std::string, std::pair<int,int> > periods;
2584
2585 ts.GetLHCPeriodBoundaries(periods);
2586
2587 TLegend* legend = new TLegend(0.15,0.82,0.90,0.92);
2588 legend->SetFillColor(0);
2589 Int_t n(0);
2590
2591
2592 for ( std::map<std::string, std::pair<int,int> >::const_iterator pit = periods.begin(); pit != periods.end(); ++pit )
2593 {
2594 std::string period = pit->first;
2595 int run1 = (pit->second).first;
2596 int run2 = (pit->second).second;
2597 int nruns(0);
2598 for ( std::set<int>::const_iterator rit = RunNumbers().begin(); rit != RunNumbers().end(); ++ rit )
2599 {
2600 if ( (*rit) >= run1 && (*rit) <= run2 )
2601 {
2602 ++nruns;
2603 }
2604 }
2605 AliInfo(Form("Period %s runs %6d-%6d ; %d actual runs",period.c_str(),run1,run2,nruns));
2606
2607 g->Fit("pol0","+Q","",run1,run2);
2608 TF1* func = static_cast<TF1*>(g->GetListOfFunctions()->Last());
2609 if (func)
2610 {
2611 func->SetLineColor(2+n);
2612 legend->AddEntry(func,Form("%s %5.2f #pm %5.2f %s (rel. error %5.2f %%)",period.c_str(),func->GetParameter(0),func->GetParError(0),
2613 (asRejection ? "%":""),100*func->GetParError(0)/func->GetParameter(0)));
2614 ++n;
2615 }
2616 }
2617
2618 legend->SetNColumns(3);
2619
2620 Double_t mean = TMath::Mean(g->GetN(),g->GetY());
2621 Double_t rms = TMath::RMS(g->GetN(),g->GetY());
2622
2623 legend->AddEntry("",Form("Mean %5.2f RMS %5.2f (%5.2f %%)",mean,rms,(mean) ? 100.0*rms/mean : 0.0),"");
2624
2625 legend->Draw();
2626
2627 return g;
2628}
2629
2630//_____________________________________________________________________________
2631void AliAnalysisMuMu::ShowList(const char* title, const TString& list, const char separator) const
2632{
2633 /// Show a list of strings
2634 TObjArray* parts = list.Tokenize(separator);
2635
2636 std::cout << title << " (" << parts->GetEntries() << ") : " << std::endl;
2637
2638 TIter next(parts);
2639 TObjString* str;
2640
2641 while ( ( str = static_cast<TObjString*>(next()) ) )
2642 {
2643 std::cout << " " << str->String().Data() << std::endl;
1afce1ce 2644 }
2645
cb690a12 2646 if ( parts->GetEntries()==0) std::cout << endl;
2647
2648 delete parts;
2754fd07 2649}
2650
2f331ac9 2651//_____________________________________________________________________________
cb690a12 2652void AliAnalysisMuMu::Print(Option_t* opt) const
2f331ac9 2653{
cb690a12 2654 /// printout
2655 std::cout << "Reading from file : " << fFilename.Data() << std::endl;
2f331ac9 2656
cb690a12 2657 TString copt(opt);
2658
2659 if (IsSimulation() || SIM() )
2f331ac9 2660 {
cb690a12 2661 copt += "SIM";
2f331ac9 2662 }
cb690a12 2663
2664 if ( !IsSimulation() )
2f331ac9 2665 {
cb690a12 2666 copt += "REAL";
2f331ac9 2667 }
2668
cb690a12 2669 Config()->Print(copt.Data());
81190958 2670
cb690a12 2671 if ( RunNumbers().size() > 1 )
81190958 2672 {
cb690a12 2673 std::cout << RunNumbers().size() << " runs";
2674 }
2675 else
2676 {
2677 std::cout << RunNumbers().size() << " run";
81190958 2678 }
2679
cb690a12 2680 if ( fCorrectionPerRun )
2681 {
2682 std::cout << " with correction factors";
2683 }
1a37b253 2684 std::cout << std::endl;
cb690a12 2685 Int_t i(0);
2686 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
81190958 2687 {
cb690a12 2688 std::cout << (*it);
2689 if ( fCorrectionPerRun )
81190958 2690 {
cb690a12 2691 std::cout << Form("(%e)",fCorrectionPerRun->GetY()[i]);
81190958 2692 }
cb690a12 2693 std::cout << ",";
2694 ++i;
81190958 2695 }
cb690a12 2696 std::cout << std::endl;
81190958 2697
cb690a12 2698 TString sopt(opt);
2699 sopt.ToUpper();
81190958 2700
cb690a12 2701 if ( sopt.Contains("BIN") && BIN() )
81190958 2702 {
cb690a12 2703 std::cout << "Binning : " << std::endl;
2704 TString topt(sopt);
2705 topt.ReplaceAll("BIN","");
2706 BIN()->Print(topt.Data());
2707 }
2708 if ( sopt.Contains("MC") && OC() )
2709 {
2710 TString topt(sopt);
2711 topt.ReplaceAll("MC","");
2712 OC()->Print(topt.Data());
2713 }
2714 if ( sopt.Contains("CC") && CC() )
2715 {
2716 CC()->Print("trigger/event");
81190958 2717 }
81190958 2718
cb690a12 2719 if ( sopt.Contains("SIZE") )
81190958 2720 {
cb690a12 2721 TFile* f = ReOpen(fFilename,"READ");
2722 TIter next(f->GetListOfKeys());
2723 TKey* key;
81190958 2724
cb690a12 2725 while ( ( key = static_cast<TKey*>(next()) ) )
81190958 2726 {
cb690a12 2727 std::cout << key->GetName() << " " << key->GetNbytes() << " " << key->GetObjlen() << std::endl;
81190958 2728 }
cb690a12 2729 }
2730}
81190958 2731
cb690a12 2732//_____________________________________________________________________________
2733TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
2734{
2735 /// Tries to reopen the file with a new mode
2736
2737 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(filename));
2738
2739 if (f)
2740 {
2741 delete f;
81190958 2742 }
2743
cb690a12 2744 f = TFile::Open(filename,mode);
81190958 2745
cb690a12 2746 if ( !f || !f->IsOpen() )
2747 {
2748 AliError(Form("Cannot open file %s in mode %s",filename,mode));
2749 return 0x0;
2750 }
2751
2752 return f;
2753}
81190958 2754
cb690a12 2755//_____________________________________________________________________________
2756void AliAnalysisMuMu::SetCentralitySelectionList(const char* centralitySelectionList)
2757{
2758 /// Set centralities to be used during fitting
2759 /// centralitySelectionList is a regular expression.
81190958 2760
cb690a12 2761 TObjArray* centralities = BIN()->CreateBinObjArray("centrality");
2762 TIter next(centralities);
2763 AliAnalysisMuMuBinning::Range* r;
2764
2765 TString csl;
2766
2767 TPRegexp re(centralitySelectionList);
2768
2769 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
81190958 2770 {
cb690a12 2771 AliDebug(1,Form("r=%s",r->AsString().Data()));
2772 if ( re.MatchB(r->AsString()) )
81190958 2773 {
cb690a12 2774 csl += r->AsString();
2775 csl += ",";
81190958 2776 }
81190958 2777 }
2778
cb690a12 2779 Config()->SetList(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation(),csl);
81190958 2780
cb690a12 2781 delete centralities;
81190958 2782}
2783
2784//_____________________________________________________________________________
cb690a12 2785Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
81190958 2786{
cb690a12 2787 /// Sets the graph used to correct values per run
2788 delete fCorrectionPerRun;
2789 fCorrectionPerRun=0x0;
81190958 2790
cb690a12 2791 // check that corr has the same runs as we do
81190958 2792
cb690a12 2793 Int_t i(0);
81190958 2794
cb690a12 2795 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
81190958 2796 {
cb690a12 2797 Int_t corrRun = TMath::Nint(corr.GetX()[i]);
2798 if (corrRun != *it)
81190958 2799 {
cb690a12 2800 AliError(Form("%d-th run mistmatch %d vs %d",i,corrRun,*it));
81190958 2801
cb690a12 2802 return kFALSE;
81190958 2803 }
cb690a12 2804 ++i;
2805 }
2806
2807 fCorrectionPerRun = new TGraphErrors(corr.GetN());
81190958 2808
cb690a12 2809 TFormula* tformula(0x0);
2810 if ( strlen(formula) > 0 )
2811 {
2812 tformula = new TFormula("SetCorrectionPerRunFormula",formula);
81190958 2813 }
2814
cb690a12 2815 i = 0;
81190958 2816
cb690a12 2817 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2818 {
2819 Double_t y = corr.GetY()[i];
2820
2821 if ( tformula )
2822 {
2823 y = tformula->Eval(y);
2824 }
2825 fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
2826 ++i;
2827 }
81190958 2828
cb690a12 2829 delete formula;
81190958 2830
cb690a12 2831 return kTRUE;
2832}
81190958 2833
cb690a12 2834//_____________________________________________________________________________
2835void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuJpsiResult& r)
2836{
2837 /// Set the "NofInput" variable(s) of one result
81190958 2838
cb690a12 2839 TString hname(Form("MinvUS%s",r.Bin().AsString().Data()));
2840
2841 TH1* hinput = fMergeableCollection->Histo(Form("/%s/ALL/ANY/V0A/INYRANGE",AliAnalysisMuMuBase::MCInputPrefix()),hname.Data());
2842
2843 if (!hinput)
2844 {
2845 AliError(Form("Got a simulation file where I did not find histogram /%s/ALL/EVERYTHING/ANY/INYRANGE/%s",AliAnalysisMuMuBase::MCInputPrefix(),hname.Data()));
81190958 2846
cb690a12 2847 }
2848 else
2849 {
2850 r.SetNofInputParticles(*hinput);
2851 }
81190958 2852}
2853
2f331ac9 2854//_____________________________________________________________________________
cb690a12 2855AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
2f331ac9 2856{
cb690a12 2857 /// Shortcut method to get to a spectra
2858 if (!OC()) return 0x0;
2f331ac9 2859
cb690a12 2860 return static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(fullpath));
2f331ac9 2861}
2862
a58729a5 2863//_____________________________________________________________________________
cb690a12 2864void AliAnalysisMuMu::SelectRunByTrigger(const char* triggerList)
a58729a5 2865{
cb690a12 2866 if (!fMergeableCollection || !fCounterCollection) return;
2867
2868 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2869 TIter nextRun(runs);
2870
2871 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2872 TIter nextTrigger(triggers);
2873
2874 TObjString* srun;
2875 TObjString* strigger;
a58729a5 2876
cb690a12 2877 TString striggerList(triggerList);
a58729a5 2878
cb690a12 2879 TList* runList = new TList();
a58729a5 2880
cb690a12 2881 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2882 {
2883
2884 nextTrigger.Reset();
2885
2886 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2887 {
2888 if ( !striggerList.Contains(strigger->String().Data()) )
2889 {
2890 continue;
2891 }
2892
2893 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2894 strigger->String().Data(),"PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00",srun->String().Atoi())));
2895 if ( n > 0 ) runList->Add(srun);
2896
2897 }
2898 }
2899 runList->Sort();
2900 TIter nextRunOK(runList);
2901 while ( ( srun = static_cast<TObjString*>(nextRunOK()) ) )
2902 {
2903 std::cout << srun->String().Atoi() << std::endl;
2904 }
2905 delete runList;
a58729a5 2906
cb690a12 2907 delete triggers;
2908 delete runs;
2909}
a58729a5 2910//_____________________________________________________________________________
cb690a12 2911void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
2912 Bool_t compact,
2913 Bool_t orderByTriggerCount)
a58729a5 2914{
cb690a12 2915 // Give the fraction of triggers (in triggerList) relative
2916 // to what is expected in the scalers
a58729a5 2917
cb690a12 2918 TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
a58729a5 2919
cb690a12 2920 AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
a58729a5 2921
cb690a12 2922 AliLog::SetGlobalLogLevel(AliLog::kFatal);
a58729a5 2923
cb690a12 2924 if (!fMergeableCollection || !fCounterCollection) return;
a58729a5 2925
81190958 2926 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
cb690a12 2927 TIter nextRun(runs);
a58729a5 2928
cb690a12 2929 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2930 TIter nextTrigger(triggers);
a58729a5 2931
cb690a12 2932 TObjString* srun;
2933 TObjString* strigger;
1afce1ce 2934
cb690a12 2935 TString striggerList(triggerList);
a58729a5 2936
cb690a12 2937 ULong64_t total(0);
2938 ULong64_t totalExpected(0);
2939 TString msg;
2940 std::multimap<ULong64_t,std::string> messages;
2941
2942 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
a58729a5 2943 {
cb690a12 2944 msg.Form("RUN %09d ",srun->String().Atoi());
a58729a5 2945
cb690a12 2946 if (!compact)
2947 {
2948 msg += "\n";
2949 }
2950
2951 ULong64_t nmax(0);
2952
2953 nextTrigger.Reset();
a58729a5 2954
cb690a12 2955 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
a58729a5 2956 {
cb690a12 2957 if ( !striggerList.Contains(strigger->String().Data()) )
a58729a5 2958 {
a58729a5 2959 continue;
2960 }
cb690a12 2961
2962 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2963 strigger->String().Data(),"ALL",srun->String().Atoi())));
a58729a5 2964
cb690a12 2965 ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
1afce1ce 2966
1afce1ce 2967
cb690a12 2968 nmax = TMath::Max(n,nmax);
2969
2970 total += n;
2971 totalExpected += expected;
2972
2973 msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
2974 (n>expected ? "!" : " "));
2975
2976 if ( expected > 0 ) {
2977 msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
1afce1ce 2978 }
cb690a12 2979
2980 if (!compact)
1afce1ce 2981 {
cb690a12 2982 msg += "\n";
1afce1ce 2983 }
a58729a5 2984 }
cb690a12 2985 if (nmax>0)
a58729a5 2986 {
cb690a12 2987 if (!orderByTriggerCount)
1afce1ce 2988 {
cb690a12 2989 std::cout << msg.Data() << std::endl;
2990 }
2991 else
2992 {
2993 messages.insert(std::make_pair(nmax,static_cast<std::string>(msg.Data())));
1afce1ce 2994 }
a58729a5 2995 }
cb690a12 2996 }
a58729a5 2997
cb690a12 2998 std::multimap<ULong64_t,std::string>::const_reverse_iterator it;
a58729a5 2999
cb690a12 3000 ULong64_t current(0);
3001 Int_t n(0);
a58729a5 3002
cb690a12 3003 for ( it = messages.rbegin(); it != messages.rend(); ++it )
3004 {
3005 ++n;
3006 current += it->first;
3007 Double_t percent = ( total > 0.0 ? current*100.0/total : 0.0);
3008 std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",percent,n ) << std::endl;
3009 }
3010
3011 std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
3012 total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
3013
3014
3015
3016 AliLog::SetGlobalLogLevel(oldLevel);
3017 delete triggers;
3018 delete runs;
a58729a5 3019}
3020
1afce1ce 3021//_____________________________________________________________________________
cb690a12 3022void AliAnalysisMuMu::UnsetCorrectionPerRun()
1afce1ce 3023{
cb690a12 3024 // drop the correction factors
3025 delete fCorrectionPerRun;
3026 fCorrectionPerRun=0x0;
3027}
1afce1ce 3028
cb690a12 3029//_____________________________________________________________________________
3030void AliAnalysisMuMu::Update()
3031{
3032 /// update the current file with memory
3033
3034 if (!CC() || !OC()) return;
1afce1ce 3035
cb690a12 3036 ReOpen(fFilename,"UPDATE");
3037
3038 if (OC())
1afce1ce 3039 {
cb690a12 3040 OC()->Write("OC",TObject::kSingleKey|TObject::kOverwrite);
1afce1ce 3041 }
cb690a12 3042
3043 ReOpen(fFilename,"READ");
1afce1ce 3044
cb690a12 3045 GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
3046}
3047
3048//_____________________________________________________________________________
3049Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
3050{
3051 /// Upgrade a file
3052 AliAnalysisMuMu m(filename);
1afce1ce 3053
cb690a12 3054 return m.Upgrade();
1afce1ce 3055}
3056
2f331ac9 3057//_____________________________________________________________________________
cb690a12 3058Bool_t AliAnalysisMuMu::Upgrade()
2f331ac9 3059{
cb690a12 3060 /// Upgrade the current file
3061 /// - from single list to one key per object, if needed
3062 /// - from histogramCollection to mergeableCollection, if needed
3063
2f331ac9 3064
cb690a12 3065 AliWarning("Out of date method");
2f331ac9 3066
cb690a12 3067 TFile* f = ReOpen(fFilename,"UPDATE");
2f331ac9 3068
cb690a12 3069 TList* list = static_cast<TList*>(f->Get("chist"));
3070
3071 if (list)
3072 {
3073 // really old file where everything was in a single list
3074
3075 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(list->At(0));
3076 AliCounterCollection* cc = static_cast<AliCounterCollection*>(list->At(1));
3077
3078 AliMergeableCollection* mc = hc->Convert();
3079
3080 f->cd();
3081
3082 mc->Write("MC",TObject::kSingleKey);
3083 cc->Write("CC",TObject::kSingleKey);
3084
3085 f->Delete("chist;*");
3086
3087 f->Write();
3088
3089 }
3090 else
2f331ac9 3091 {
cb690a12 3092 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(f->Get("HC"));
3093
3094 if ( hc )
2f331ac9 3095 {
cb690a12 3096 // old file with histogram collection instead of mergeable collection
3097
3098 AliMergeableCollection* mc = hc->Convert();
3099
3100 f->cd();
3101
3102 mc->Write("MC",TObject::kSingleKey);
3103
3104 f->Delete("HC;*");
3105
3106 f->Write();
2f331ac9 3107 }
2f331ac9 3108 }
cb690a12 3109
3110 delete f;
2f331ac9 3111
cb690a12 3112 return kTRUE;
2f331ac9 3113}
3114
1afce1ce 3115//_____________________________________________________________________________
cb690a12 3116TH2* AliAnalysisMuMu::ComputeSPDCorrection(const char* type, const char* eventSel, const char* triggerSel, Bool_t bkgReject)
1afce1ce 3117{
cb690a12 3118 TString stype(type);
3119 TString evtype(eventSel);
3120 TString trigtype(triggerSel);
1afce1ce 3121
cb690a12 3122 TH2* hNch = static_cast<TH2*>(OC()->Histo(Form("/MCINPUT/%s/%s/V0A/NchVsZVertexVsEta",evtype.Data(),
3123 trigtype.Data()))); // Input Nch // //"/INPUT/QASPDZPSALL/NchVSEtaVSZVertMC"
3124 if ( !hNch )
1afce1ce 3125 {
cb690a12 3126 AliError("No Nch histo found");
3127 return 0x0;
1afce1ce 3128 }
cb690a12 3129 TH2* hNtr = static_cast<TH2*>(OC()->Histo(Form("/%s/%s/V0A/TrackletsVsZVertexVsEta",evtype.Data(),
3130 trigtype.Data()))); // Reco tracklets // //"/RECO/QASPDZPSALL/MB1/NtrVSEtaVSZVertMC"
3131 if ( !hNtr )
1afce1ce 3132 {
cb690a12 3133 AliError("No tracklets histo found");
3134 return 0x0;
1afce1ce 3135 }
2f331ac9 3136
cb690a12 3137 TH2* hNtrBkg = static_cast<TH2*>(OC()->Histo(Form("/MCINPUT/%s/%s/V0A/NBkgTrackletsVsZVertexVsEta",evtype.Data(),
3138 trigtype.Data()))); // Reco tracklets // //"/RECO/QASPDZPSALL/MB1/NtrVSEtaVSZVertMC"
3139 if ( !hNtrBkg )
3140 {
3141 AliError("No background tracklets histo found");
3142 return 0x0;
3143 }
3144
2f331ac9 3145
cb690a12 3146 TH2D* hSPDCorr = static_cast<TH2D*>(hNtr->Clone("SPDCorr"));
3147 TString title("");\
3148 if ( stype.Contains("oneOverAccEff")) hSPDCorr->SetTitle("SPD 1/AccxEff correction");
3149 else if ( stype.Contains("AccEffOnly")) hSPDCorr->SetTitle("SPD AccxEff correction");
3150 else if ( stype.Contains("statOneOverAccEff")) hSPDCorr->SetTitle("SPD 1/AccxEff correction stat. unc.");
2f331ac9 3151
cb690a12 3152 for (Int_t i = 1 ; i < hNch->GetNbinsX() ; i++)
2f331ac9 3153 {
cb690a12 3154 for (Int_t j = 1 ; j < hNch->GetNbinsY() ; j++)
2f331ac9 3155 {
cb690a12 3156 Int_t n = hNch->GetBin(i,j);
3157 Double_t nch = hNch->GetBinContent(n);
3158 Double_t ntr = hNtr->GetBinContent(n);
3159 Double_t nBkgtr(0.);
3160 if ( bkgReject ) nBkgtr = hNtrBkg->GetBinContent(n);
3161
3162 Double_t corr(0.),corrErr(0.);
3163 if ( nch != 0. )
2f331ac9 3164 {
cb690a12 3165 corr = (ntr - nBkgtr)/nch;
3166 corrErr = TMath::Max( 1./nch,TMath::Sqrt( corr*(1.-corr)/nch ) );
3167 }
3168
3169 if ( stype.Contains("oneOverAccEff"))
3170 {
3171 if ( corr > 0. )
2f331ac9 3172 {
cb690a12 3173 hSPDCorr->SetBinContent(n,1./corr);
3174 hSPDCorr->SetBinError(n,corrErr/TMath::Power(corr,2.));
2f331ac9 3175 }
3176 else
3177 {
cb690a12 3178 hSPDCorr->SetBinContent(n,0.);
3179 hSPDCorr->SetBinError(n,1.);
3180 }
3181
3182 }
3183 else if ( stype.Contains("AccEffOnly"))
3184 {
3185 hSPDCorr->SetBinContent(n,corr);
3186 hSPDCorr->SetBinError(n,corrErr);
3187 }
3188 else if ( stype.Contains("statOneOverAccEff"))
3189 {
3190 if ( corr != 0. )
3191 {
3192 hSPDCorr->SetBinContent(n,(corrErr/TMath::Power(corr,2.)*100)/(1./corr));
3193 }
3194 else
3195 {
3196 hSPDCorr->SetBinContent(n,-1);
2f331ac9 3197 }
cb690a12 3198
2f331ac9 3199 }
3200 }
3201 }
3202
cb690a12 3203 return hSPDCorr;
2f331ac9 3204}
3205
3206//_____________________________________________________________________________
cb690a12 3207void AliAnalysisMuMu::ComputeFnorm()
2f331ac9 3208{
cb690a12 3209 /// Compute the CMUL to CINT ratio(s)
a58729a5 3210
cb690a12 3211 if (!CC()) return;
a58729a5 3212
cb690a12 3213 OC()->Prune("/FNORM");
a58729a5 3214
cb690a12 3215 AliAnalysisMuMuFnorm computer(*(CC()),AliAnalysisMuMuFnorm::kMUL,Config()->OCDBPath(),Config()->CompactGraphs());
a58729a5 3216
cb690a12 3217 computer.ComputeFnorm();
a58729a5 3218
cb690a12 3219 AliMergeableCollection* fnorm = computer.DetachMC();
a58729a5 3220
cb690a12 3221 OC()->Attach(fnorm,"/FNORM/");
81190958 3222
cb690a12 3223 Update();
a58729a5 3224}
3225
3226//_____________________________________________________________________________
cb690a12 3227TH1* AliAnalysisMuMu::ComputeDiffFnormFromHistos(const char* what,const char* quantity,const char* flavour,Bool_t printout)
a58729a5 3228{
cb690a12 3229 /// Compute the CMUL to CINT ratio(s) from the histos stored in the OC()
3230 //FIXME: This is just a patch...
3231
3232 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3233 if ( !binning )
a58729a5 3234 {
cb690a12 3235 AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3236 return 0x0;
a58729a5 3237 }
cb690a12 3238 TObjArray* dNchdEtas = binning->CreateBinObjArray();
a58729a5 3239
cb690a12 3240 Double_t* binArray = binning->CreateBinArray();
a58729a5 3241
cb690a12 3242 TIter next(dNchdEtas);
3243 AliAnalysisMuMuBinning::Range* r;
a58729a5 3244
cb690a12 3245 Double_t FNorm(0.);
3246 Double_t FNormError(0.);
a58729a5 3247
cb690a12 3248 TH1* hFNorm = new TH1F("hFNorm","'Global' normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",dNchdEtas->GetEntries(),binArray);
a58729a5 3249
cb690a12 3250 Int_t bin(0);
3251 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
3252 {
a58729a5 3253
cb690a12 3254 TH1* hCMSL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMSL7-B-NOPF-MUON/V0A/%s",
3255 Form("EventsIn%s",r->AsString().Data())));
3256 if ( !hCMSL )
a58729a5 3257 {
cb690a12 3258 AliError(Form("No event histo in bin %s found for CMSL7-B-NOPF-MUON",r->AsString().Data()));
3259 delete binning;
3260 delete dNchdEtas;
3261 delete binArray;
3262 delete hFNorm;
3263 return 0x0;
a58729a5 3264 }
cb690a12 3265
3266 TH1* hCMSLandOMUL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMSL7-B-NOPF-MUON&0MUL/V0A/%s",
3267 Form("EventsIn%s",r->AsString().Data())));
3268 if ( !hCMSLandOMUL )
a58729a5 3269 {
cb690a12 3270 AliError(Form("No event histo in bin %s found for CMSL7-B-NOPF-MUON & 0MUL",r->AsString().Data()));
3271 delete binning;
3272 delete dNchdEtas;
3273 delete binArray;
3274 delete hFNorm;
3275 return 0x0;
a58729a5 3276 }
2f331ac9 3277
cb690a12 3278 TH1* hCINT = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CINT7-B-NOPF-ALLNOTRD/V0A/%s",
3279 Form("EventsIn%s",r->AsString().Data())));
3280 if ( !hCINT )
a58729a5 3281 {
cb690a12 3282 AliError(Form("No event histo in bin %s found for CINT7-B-NOPF-ALLNOTRD",r->AsString().Data()));
3283 delete binning;
3284 delete dNchdEtas;
3285 delete binArray;
3286 delete hFNorm;
3287 return 0x0;
a58729a5 3288 }
81190958 3289
cb690a12 3290 TH1* hCINTandOMSL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CINT7-B-NOPF-ALLNOTRD&0MSL/V0A/%s",
3291 Form("EventsIn%s",r->AsString().Data())));
3292 if ( !hCINTandOMSL )
2f331ac9 3293 {
cb690a12 3294 AliError(Form("No event histo in bin %s found for CINT7-B-NOPF-ALLNOTRD & 0MSL",r->AsString().Data()));
3295 delete binning;
3296 delete dNchdEtas;
3297 delete binArray;
3298 delete hFNorm;
3299 return 0x0;
2f331ac9 3300 }
cb690a12 3301
3302 FNorm = (hCMSL->GetBinContent(1)/hCMSLandOMUL->GetBinContent(1))*(hCINT->GetBinContent(1)/hCINTandOMSL->GetBinContent(1));
3303 FNormError = ErrorPropagationAxBoverCxD(hCMSL->GetBinContent(1),hCINT->GetBinContent(1),hCMSLandOMUL->GetBinContent(1),hCINTandOMSL->GetBinContent(1));
3304
3305 if ( printout ) std::cout << r->AsString().Data() << " : " << FNorm << " +- " << FNormError << std::endl;
3306
3307 hFNorm->SetBinContent(++bin,FNorm);
3308 hFNorm->SetBinError(bin,FNormError);
3309 }
3310
3311 delete binning;
3312 delete dNchdEtas;
3313 delete[] binArray;
3314
3315 return hFNorm;
3316}
3317
3318//_____________________________________________________________________________
3319void AliAnalysisMuMu::ComputeDiffFnormFromInt(const char* triggerCluster, const char* eventSelection, AliMergeableCollection* mc, const char* what,const char* quantity,const char* flavour,Bool_t printout)
3320{
3321 TString striggerCluster(triggerCluster);
3322 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3323 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3324 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3325 else
3326 {
3327 AliError("Unknown trigger cluster");
3328 return;
3329 }
3330
3331 TString seventSelection(eventSelection);
3332
3333 TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
3334
3335 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3336 if ( !binning )
3337 {
3338 AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3339 return;
3340 }
3341
3342 TString path(Form("%s/%s/%s",
3343 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
3344 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
3345 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
3346 if ( !mc )
3347 {
3348 AliError("Error: No mergeable collection to get Nch histo");
3349 delete binning;
3350 return;
3351 }
3352 TH1* hNtr = static_cast<TH1*>(mc->Histo(Form("/%s/Nch",path.Data())));
3353 if ( !hNtr )
3354 {
3355 AliError(Form("Error: No /%s/Nch histo in mergeable collection",path.Data()));
3356 delete binning;
3357 return;
2f331ac9 3358 }
cb690a12 3359 Int_t nTrackletsCorrTot = hNtr->Integral();
2f331ac9 3360
cb690a12 3361 TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3362 Int_t nEntries = bin->GetEntries();
3363 Double_t* binArray = binning->CreateBinArray();
3364 Double_t FNormTot(0.);
3365 Double_t FNormTotError(0.);
3366
3367 TH1* hNorm = OC()->Histo(Form("%s/hFNormInt",id.Data()));
3368
3369 TH1* hFNormTot = new TH1F("hFNormVSdNchdEtaFromInt","Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",nEntries,binArray);
3370
3371 Double_t FNormGlobal = hNorm->GetBinContent(1);
3372 Double_t FNormGlobalError = hNorm->GetBinError(1);
3373
3374 if ( printout ) std::cout << "Global FNorm = " << FNormGlobal << " + - " << FNormGlobalError << std::endl;
3375
3376 TIter nextBin(bin);
3377 AliAnalysisMuMuBinning::Range* r;
3378 Int_t i(1);
3379 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
3380 {
3381 Double_t nTrklsCorrBin = hNtr->Integral(r->Xmin(),r->Xmax());
3382 Double_t nTrklsCorrBinFrac = nTrklsCorrBin / nTrackletsCorrTot;
3383 Double_t nTrklsCorrBinFracError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nTrklsCorrBin) / nTrackletsCorrTot,2.) +
3384 TMath::Power(TMath::Sqrt(nTrackletsCorrTot)*nTrklsCorrBin / TMath::Power(nTrackletsCorrTot,2.) ,2.) );
3385
3386 FNormTot = FNormGlobal*nTrklsCorrBinFrac;
3387 FNormTotError = TMath::Sqrt( TMath::Power(FNormGlobalError*nTrklsCorrBinFrac,2.) + TMath::Power(FNormGlobal*nTrklsCorrBinFracError,2.) );
1afce1ce 3388
cb690a12 3389 hFNormTot->SetBinContent(i,FNormTot);
3390 hFNormTot->SetBinError(i,FNormTotError);
3391 i++;
3392
3393 if ( printout ) std::cout << "Bin: " << r->AsString().Data() << " ; " << " FNorm = " << FNormTot << " +- " << FNormTotError << std::endl;
3394 }
3395
3396 TH1* o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
3397
3398 if (o)
2f331ac9 3399 {
cb690a12 3400 AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
3401 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
a58729a5 3402 }
1afce1ce 3403
cb690a12 3404 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
1afce1ce 3405
cb690a12 3406 if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
3407 else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
a58729a5 3408
cb690a12 3409 delete binning;
3410 delete bin;
3411 delete[] binArray;
a58729a5 3412}
3413
81190958 3414//_____________________________________________________________________________
cb690a12 3415void AliAnalysisMuMu::ComputeDiffFnormFromCounters(const char* triggerCluster, const char* eventSelection, const char* filePileUpCorr, const char* what,const char* quantity,const char* flavour,Bool_t printout)
81190958 3416{
cb690a12 3417 /// Compute the CMUL to CINT ratio(s) in 2 steps from the CC(), in bins
81190958 3418
cb690a12 3419 TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3420 if ( colType.Contains("-B-") ) colType = "B";
3421 else if ( colType.Contains("-S-") ) colType = "S";
3422 else
3423 {
3424 AliError("Unknown collision type");
3425 return;
3426 }
81190958 3427
cb690a12 3428 TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3429 if ( triggerType.Contains("7-") ) triggerType = "7";
3430 else if ( triggerType.Contains("8-") ) triggerType = "8";
3431 else
3432 {
3433 AliError("Unknown trigger type");
3434 return;
3435 }
81190958 3436
cb690a12 3437 TString striggerCluster(triggerCluster);
3438 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3439 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3440 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3441 else
3442 {
3443 AliError("Unknown trigger cluster");
3444 return;
3445 }
81190958 3446
cb690a12 3447 std::cout << striggerCluster.Data() << std::endl;
81190958 3448
cb690a12 3449 Bool_t corrPU(kFALSE);
3450 TObjArray* pUCorr = new TObjArray();
3451 if ( strlen(filePileUpCorr) > 0 )
81190958 3452 {
cb690a12 3453 // std::cout << "Extracting Pile-Up correction factors from " << filePileUpCorr << std::endl;
3454 char line[1024];
3455 ifstream in(filePileUpCorr);
81190958 3456
cb690a12 3457 while ( in.getline(line,1024,'\n'))
81190958 3458 {
cb690a12 3459 TString lrun(line);
3460 TString lvalue(line);
3461
3462 lrun.Remove(0,4);
3463 lrun.Remove(6,67);
3464
3465 lvalue.Remove(0,57);//71
3466
3467 // std::cout << "RUN: " << lrun.Data() << " PUFactor = " << lvalue.Data() << std::endl;
3468
3469 pUCorr->Add(new TParameter<Double_t>(lrun.Data(),lvalue.Atof()));
81190958 3470 }
cb690a12 3471 corrPU = kTRUE;
81190958 3472 }
81190958 3473
cb690a12 3474 TString seventSelection(eventSelection);
3475 TString sruns = CC()->GetKeyWords("run");
3476 TObjArray* runs = sruns.Tokenize(",");
3477 Double_t NofRuns = runs->GetEntries();
3478
3479 TIter nextRun(runs);
3480 TObjString* s;
3481
3482 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3483 if ( !binning )
81190958 3484 {
cb690a12 3485 AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3486 return;
3487 }
3488 TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3489 Double_t* binArray = binning->CreateBinArray();
3490 Int_t nEntries = bin->GetEntries();
3491
3492 TH1* h;
3493 TH1* hNofEqMB = new TH1F("hNofEqMBVSdNchdEta","Equivalent MB events per CMUL for vs dN_{ch}/d#eta",bin->GetEntries(),binArray);
3494 TH1* hFNormTot = new TH1F("hFNormVSdNchdEta","Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",bin->GetEntries(),binArray);
81190958 3495
cb690a12 3496 Double_t* FNormTot = new Double_t[nEntries];
3497 Double_t* FNormTotError = new Double_t[nEntries];
81190958 3498
cb690a12 3499 TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3500
3501 TList* lRun2Reject = new TList();
3502 lRun2Reject->SetOwner(kTRUE);
81190958 3503
cb690a12 3504 Int_t i(0); //dNchdEta bin number
3505 TObjArray* aCluster = striggerCluster.Tokenize("-");
3506 TIter nextCluster(aCluster);
3507 TObjString* striggerClusterS;
3508
3509 TIter nextBin(bin);
3510 AliAnalysisMuMuBinning::Range* r;
3511 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
3512 {
3513 FNormTot[i] = 0;
3514 FNormTotError[i] = 0;
3515 if ( printout )
81190958 3516 {
cb690a12 3517 std::cout << "______________________________" << std::endl;
3518 std::cout << "Bin: " << r->AsString().Data() << std::endl;
81190958 3519 }
81190958 3520
cb690a12 3521 h = new TH1F(Form("hFNormVSrun_%s",r->AsString().Data()),Form("Normalization factor vs run for %s ;run;FNorm",r->AsString().Data()),NofRuns,1,NofRuns);
3522 //Set the run labels
81190958 3523
cb690a12 3524 Double_t nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3525 seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data())); //Nof CMUL7/8 events in Bin summed over runs
3526 //HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
81190958 3527
cb690a12 3528 Int_t j(1); //Run label index
3529 nextRun.Reset();
3530 while ( ( s = static_cast<TObjString*>(nextRun())) ) //Run loop
81190958 3531 {
cb690a12 3532 Double_t nCMSL(0.),nCMSLandOMUL(0.);
3533 nextCluster.Reset();
3534 while ( (striggerClusterS = static_cast<TObjString*>(nextCluster())) && nCMSL == 0. )
81190958 3535 {
cb690a12 3536 nCMSL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s/centrality:V0A/run:%s/bin:%s",
3537 seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName(),r->AsString().Data()));
81190958 3538
cb690a12 3539 nCMSLandOMUL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s&0MUL/centrality:V0A/run:%s/bin:%s",
3540 seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName(),r->AsString().Data()));
81190958 3541 }
cb690a12 3542 Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s/bin:%s",
3543 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3544
3545 Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/run:%s/bin:%s",
3546 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3547
3548 Double_t nCINTandOMSL = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD&0MSL/centrality:V0A/run:%s/bin:%s",
3549 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3550
3551 Double_t FNorm(0.);
3552 Double_t FNormError(0.);
3553 Double_t FNormError2(0.);
3554 Double_t pUfactor = 1.;
3555 if ( nCMSLandOMUL != 0. && nCINTandOMSL !=0. && nCMSL != 0. && nCINT !=0. )
81190958 3556 {
cb690a12 3557 if (corrPU)
3558 {
3559 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(pUCorr->FindObject(s->GetName()));
3560 if ( p ) pUfactor = p->GetVal();
3561 else
3562 {
3563 AliError(Form("Run %s not found in pile-up correction list",s->GetName()));
3564 }
3565 }
81190958 3566
cb690a12 3567 FNorm = (nCMSL*nCINT)*pUfactor/(nCMSLandOMUL*nCINTandOMSL);
3568 FNormError = ErrorPropagationAxBoverCxD(nCMSL,nCINT,nCMSLandOMUL,nCINTandOMSL)*pUfactor;
3569 FNormError2 = AliAnalysisMuMuResult::ErrorABCD(nCMSL, TMath::Sqrt(nCMSL), nCINT, TMath::Sqrt(nCINT), nCMSLandOMUL, TMath::Sqrt(nCMSLandOMUL), nCINTandOMSL, TMath::Sqrt(nCINTandOMSL));
81190958 3570 }
3571 else
3572 {
cb690a12 3573 if ( nCINT == 0 ) std::cout << " Warning: Run " << s->GetName() << " has no MB trigger in this bin" << std::endl;
3574
3575 lRun2Reject->Add(new TObjString(s->GetName()));
3576 if ( printout ) std::cout << "Run " << s->GetName() << " not used for FNorm cause lack of stats" << std::endl;
3577 continue;
81190958 3578 }
cb690a12 3579 FNormTot[i] += FNorm*nCMUL; // This is the sum of equivalent Nof MB per CMUL run by run. NOTE: This sum is NOT always the total equivalent Nof MB per CMUL because in pp 2012 if just one cluster is used at a time this sum is not the sum for all runs
3580 FNormTotError[i] += TMath::Power(nCMUL*FNormError,2.) + TMath::Power(FNorm*TMath::Sqrt(nCMUL),2.);
81190958 3581
cb690a12 3582 if ( printout ) std::cout << "Run " << s->GetName() << " FNorm = " << FNorm << " +- " << FNormError << " (" << FNormError2 << ")" << " ; PUFactor =" << pUfactor << " ; " << "Nof CMUL = " << nCMUL << std::endl;
81190958 3583
cb690a12 3584 h->GetXaxis()->SetBinLabel(j,s->GetName());
3585 h->SetBinContent(j,FNorm);
3586 h->SetBinError(j++,FNormError);
81190958 3587
cb690a12 3588 }
3589
3590// std::cout << "NofCMUL in " << i << " = " << nCMULBin << std::endl;
3591
3592 TIter nextRejectRun(lRun2Reject);
3593 TObjString* run2Rej;
3594 Double_t nCMULBinRej(0.);
3595 while ( (run2Rej = static_cast<TObjString*>(nextRejectRun())) )
81190958 3596 {
cb690a12 3597 nCMULBinRej += CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s/run:%s",
3598 seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data(),
3599 run2Rej->GetName())); //Sum of CMUL7 events from rejected runs
81190958 3600 }
cb690a12 3601
3602 nCMULBin = nCMULBin - nCMULBinRej;
3603 lRun2Reject->Clear();
3604
3605 FNormTotError[i] = TMath::Sqrt(TMath::Power(TMath::Sqrt(FNormTotError[i])/nCMULBin,2.) + TMath::Power(FNormTot[i]*TMath::Sqrt(nCMULBin)/TMath::Power(nCMULBin,2.),2.));
3606 FNormTot[i] = FNormTot[i]/nCMULBin;
3607
3608 std::cout << "Mean FNorm in Bin = " << FNormTot[i] << " +- " << FNormTotError[i] <<std::endl;
3609
3610 hFNormTot->SetBinContent(i+1,FNormTot[i]);
3611 hFNormTot->SetBinError(i+1,FNormTotError[i]);
3612
3613 //____
3614 nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3615 seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data())); //Nof CMUL7/8 events in Bin summed over runs
3616
3617 Double_t nofEqMB = FNormTot[i]*nCMULBin;
3618 Double_t nofEqMBError = TMath::Sqrt( TMath::Power(FNormTotError[i]*nCMULBin,2.) + TMath::Power(FNormTot[i]*TMath::Sqrt(nCMULBin),2.) );
3619
3620 std::cout << "EqMB in Bin = " << nofEqMB << " +- " << nofEqMBError << " ; nCMUL = " << nCMULBin << std::endl;
3621
3622 hNofEqMB->SetBinContent(i+1,nofEqMB);
3623 hNofEqMB->SetBinError(i+1,nofEqMBError);
3624 //____
3625
3626 TH1* o = fMergeableCollection->Histo(id.Data(),h->GetName());
3627
3628 if (o)
3629 {
3630 AliWarning(Form("Replacing %s/%s",id.Data(),h->GetName()));
3631 fMergeableCollection->Remove(Form("%s/%s",id.Data(),h->GetName()));
3632 }
3633
3634 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),h);
3635
3636 if ( adoptOK ) std::cout << "+++FNorm histo " << h->GetName() << " adopted" << std::endl;
3637 else AliError(Form("Could not adopt FNorm histo %s",h->GetName()));
3638
3639 i++;
3640 lRun2Reject->Clear();
81190958 3641 }
3642
cb690a12 3643 TH1* o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
a58729a5 3644
cb690a12 3645 if (o)
a58729a5 3646 {
cb690a12 3647 AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
3648 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
a58729a5 3649 }
3650
cb690a12 3651 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
a58729a5 3652
cb690a12 3653 if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
3654 else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
a58729a5 3655
a58729a5 3656
cb690a12 3657 o = fMergeableCollection->Histo(id.Data(),hNofEqMB->GetName());
a58729a5 3658
cb690a12 3659 if (o)
3660 {
3661 AliWarning(Form("Replacing %s/%s",id.Data(),hNofEqMB->GetName()));
3662 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNofEqMB->GetName()));
3663 }
a58729a5 3664
cb690a12 3665 adoptOK = fMergeableCollection->Adopt(id.Data(),hNofEqMB);
a58729a5 3666
cb690a12 3667 if ( adoptOK ) std::cout << "+++FNorm histo " << hNofEqMB->GetName() << " adopted" << std::endl;
3668 else AliError(Form("Could not adopt FNorm histo %s",hNofEqMB->GetName()));
a58729a5 3669
cb690a12 3670 delete binning;
3671 delete runs;
3672 delete aCluster;
3673 delete bin;
3674 delete[] binArray;
3675 delete[] FNormTot;
3676 delete[] FNormTotError;
3677 delete lRun2Reject;
3678
3679 return;
3680
3681}
3682
3683//_____________________________________________________________________________
3684void AliAnalysisMuMu::ComputeDiffFnormFromGlobal(const char* triggerCluster, const char* eventSelection, const char* what,const char* quantity,
3685 const char* flavour, Bool_t printout)
3686{
3687 TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3688 if ( colType.Contains("-B-") ) colType = "B";
3689 else if ( colType.Contains("-S-") ) colType = "S";
3690 else
3691 {
3692 AliError("Unknown collision type");
3693 return;
3694 }
a58729a5 3695
cb690a12 3696 TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3697 if ( triggerType.Contains("7-") ) triggerType = "7";
3698 else if ( triggerType.Contains("8-") ) triggerType = "8";
3699 else
3700 {
3701 AliError("Unknown trigger type");
3702 return;
3703 }
a58729a5 3704
cb690a12 3705 TString striggerCluster(triggerCluster);
3706 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3707 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3708 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3709 else
a58729a5 3710 {
cb690a12 3711 AliError("Unknown trigger cluster");
3712 return;
a58729a5 3713 }
3714
cb690a12 3715 TString seventSelection(eventSelection);
3716 TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
a58729a5 3717
cb690a12 3718 TH1* hFnormGlobal = OC()->Histo(id.Data(),"hFNormVSdNchdEta");
3719 if( !hFnormGlobal)
3720 {
3721 AliError("hFNormVSdNchdEta not found");
3722 return;
3723 }
a58729a5 3724
cb690a12 3725 TH1* hFnormGlobalInt = OC()->Histo(id.Data(),"hFNormInt");
3726 if( !hFnormGlobalInt)
3727 {
3728 AliError("hFNormInt not found");
3729 return;
3730 }
3731 Double_t FNormGlobal = hFnormGlobalInt->GetBinContent(1);
3732 Double_t FNormGlobalError = hFnormGlobalInt->GetBinError(1);
a58729a5 3733
cb690a12 3734 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3735 if ( !binning )
a58729a5 3736 {
cb690a12 3737 AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3738 return;
a58729a5 3739 }
cb690a12 3740 TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
a58729a5 3741
cb690a12 3742 TH1* hFNormVSNtr = static_cast<TH1*>(hFnormGlobal->Clone());
3743 hFNormVSNtr->SetTitle("Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm");
3744 hFNormVSNtr->SetName("hFNormVSdNchdEtaFromGlobal");
a58729a5 3745
cb690a12 3746 TH1* hNMBVSNtr = static_cast<TH1*>(hFnormGlobal->Clone());
3747 hNMBVSNtr->SetTitle("Equivalent MB events per CMUL for vs dN_{ch}/d#eta;NMB");
3748 hNMBVSNtr->SetName("hNofEqMBVSdNchdEtaFromGlobal");
3749
3750
3751 Double_t nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
3752 seventSelection.Data(),triggerType.Data(),colType.Data()));
3753
3754 Double_t nCINTTot = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A",
3755 seventSelection.Data(),triggerType.Data(),colType.Data()));
3756 TIter nextBin(bin);
3757 AliAnalysisMuMuBinning::Range* r;
3758 Int_t i(1);
3759 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
a58729a5 3760 {
cb690a12 3761 Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3762 seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data()));
a58729a5 3763
cb690a12 3764 Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/bin:%s",
3765 seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data()));
a58729a5 3766
cb690a12 3767 Double_t f = nCMUL/nCMULTot;
3768 Double_t fError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nCMUL)/nCMULTot,2.) + TMath::Power(nCMUL*TMath::Sqrt(nCMULTot)/TMath::Power(nCMULTot,2.),2.) );
a58729a5 3769
cb690a12 3770 Double_t g = nCINT/nCINTTot;
3771 Double_t gError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nCINT)/nCINTTot,2.) + TMath::Power(nCINT*TMath::Sqrt(nCINTTot)/TMath::Power(nCINTTot,2.),2.) );
a58729a5 3772
cb690a12 3773 Double_t value = FNormGlobal*(g/f);
3774 Double_t error = TMath::Sqrt( TMath::Power(FNormGlobalError*(g/f),2.) + TMath::Power(FNormGlobal*(gError/f),2.) + TMath::Power(FNormGlobal*g*fError/TMath::Power(f,2.),2.) );
a58729a5 3775
cb690a12 3776 hFNormVSNtr->SetBinContent(i,value);
3777 hFNormVSNtr->SetBinError(i,error);
a58729a5 3778
cb690a12 3779 if (printout) std::cout << value << " +- " << error << " ; nCMUL = " << nCMUL << std::endl;
a58729a5 3780
cb690a12 3781 hNMBVSNtr->SetBinContent(i,value*nCMUL);
3782 hNMBVSNtr->SetBinError(i,TMath::Sqrt( TMath::Power(error*nCMUL,2.) + TMath::Power(value*TMath::Sqrt(nCMUL),2.) ));
3783
3784 i++;
2f331ac9 3785 }
2f331ac9 3786
cb690a12 3787 TH1* o = fMergeableCollection->Histo(id.Data(),hFNormVSNtr->GetName());
a58729a5 3788
cb690a12 3789 if (o)
2754fd07 3790 {
cb690a12 3791 AliWarning(Form("Replacing %s/%s",id.Data(),hFNormVSNtr->GetName()));
3792 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormVSNtr->GetName()));
2754fd07 3793 }
2f331ac9 3794
cb690a12 3795 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormVSNtr);
2f331ac9 3796
cb690a12 3797 if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormVSNtr->GetName() << " adopted" << std::endl;
3798 else AliError(Form("Could not adopt FNorm histo %s",hFNormVSNtr->GetName()));
2f331ac9 3799
cb690a12 3800 o = fMergeableCollection->Histo(id.Data(),hNMBVSNtr->GetName());
2754fd07 3801
cb690a12 3802 if (o)
2f331ac9 3803 {
cb690a12 3804 AliWarning(Form("Replacing %s/%s",id.Data(),hNMBVSNtr->GetName()));
3805 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNMBVSNtr->GetName()));
2f331ac9 3806 }
3807
cb690a12 3808 adoptOK = fMergeableCollection->Adopt(id.Data(),hNMBVSNtr);
2f331ac9 3809
cb690a12 3810 if ( adoptOK ) std::cout << "+++FNorm histo " << hNMBVSNtr->GetName() << " adopted" << std::endl;
3811 else AliError(Form("Could not adopt FNorm histo %s",hNMBVSNtr->GetName()));
2754fd07 3812
a58729a5 3813
cb690a12 3814}
a58729a5 3815
cb690a12 3816//_____________________________________________________________________________
3817void AliAnalysisMuMu::ComputeMeanFnorm(const char* triggerCluster, const char* eventSelection, const char* what,const char* quantity,const char* flavour, Bool_t printout)
3818{
3819 /// Compute the mean Fnorm and mean NMB from the offline and "rescaled global" methods
3820
3821 TString seventSelection(eventSelection);
3822 TString striggerCluster(triggerCluster);
3823 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3824 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3825 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3826 else
2f331ac9 3827 {
cb690a12 3828 AliError("Unknown trigger cluster");
3829 return;
2f331ac9 3830 }
a58729a5 3831
cb690a12 3832 TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3833 if ( triggerType.Contains("7-") ) triggerType = "7";
3834 else if ( triggerType.Contains("8-") ) triggerType = "8";
3835 else
a58729a5 3836 {
cb690a12 3837 AliError("Unknown trigger type");
3838 return;
a58729a5 3839 }
a58729a5 3840
cb690a12 3841 TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3842 if ( colType.Contains("-B-") ) colType = "B";
3843 else if ( colType.Contains("-S-") ) colType = "S";
3844 else
a58729a5 3845 {
cb690a12 3846 AliError("Unknown collision type");
3847 return;
a58729a5 3848 }
cb690a12 3849
3850 TH1* hMB = OC()->Histo(Form("/FNORM-%s/%s/V0A/hNofEqMBVSdNchdEta",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3851 if ( !hMB )
a58729a5 3852 {
cb690a12 3853 AliError("Histo hNofEqMBVSdNchdEta not found");
3854 return;
a58729a5 3855 }
cb690a12 3856
3857 TH1* hMBG = OC()->Histo(Form("/FNORM-%s/%s/V0A/hNofEqMBVSdNchdEtaFromGlobal",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3858 if ( !hMBG )
a58729a5 3859 {
cb690a12 3860 AliError("Histo hNofEqMBVSdNchdEtaFromGlobal not found");
3861 return;
a58729a5 3862 }
2f331ac9 3863
cb690a12 3864 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3865 if ( !binning )
2754fd07 3866 {
cb690a12 3867 AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3868 return;
2754fd07 3869 }
cb690a12 3870 TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3871
3872
3873 TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
2f331ac9 3874
cb690a12 3875 TH1* hMBMean = static_cast<TH1*>(hMBG->Clone());
3876 hMBMean->SetName("hNofEqMBVSdNchdEtaFromMean");
2f331ac9 3877
cb690a12 3878 TH1* hFnormMean = static_cast<TH1*>(hMBG->Clone());
3879 hFnormMean->SetName("hFNormVSdNchdEtaFromMean");
2f331ac9 3880
cb690a12 3881 for ( Int_t i = 1 ; i <= hMB->GetNbinsX() ; i++ )
2754fd07 3882 {
cb690a12 3883 Double_t Fn = hMB->GetBinContent(i);
3884 Double_t Fng = hMBG->GetBinContent(i);
2754fd07 3885
cb690a12 3886 Double_t FnE = hMB->GetBinError(i);
3887 Double_t FngE = hMBG->GetBinError(i);
2754fd07 3888
cb690a12 3889 Double_t meanBin = (Fn + Fng) / 2.;
3890 Double_t meanBinError = TMath::Sqrt( TMath::Power(FnE/2.,2.) + TMath::Power(FngE/2.,2.) );
3891 Double_t meanBinSys = TMath::Abs( meanBin - Fn );
3892
3893// Double_t meanBin = (Fn/TMath::Power(FnE,2.) + Fng/TMath::Power(FngE,2.)) / ( 1./TMath::Power(FnE,2.) + 1./TMath::Power(FngE,2.) );
3894// Double_t meanBinError = 1. / TMath::Sqrt( 1./TMath::Power(FnE,2.) + 1./TMath::Power(FngE,2.) );
3895// Double_t meanBinSys = TMath::Sqrt( TMath::Power(Fn - meanBin,2.)/TMath::Power(FnE,2.) + TMath::Power(Fng - meanBin,2.)/TMath::Power(FngE,2.) );
a58729a5 3896
a58729a5 3897
cb690a12 3898 hMBMean->SetBinContent(i,meanBin);
3899 hMBMean->SetBinError(i,meanBinError);
2754fd07 3900
cb690a12 3901 Double_t nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3902 seventSelection.Data(),triggerType.Data(),colType.Data(),
3903 static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(i-1))->AsString().Data()));
2754fd07 3904
cb690a12 3905 if (printout) std::cout << meanBinSys/nCMULBin << std::endl;
3906
3907 Double_t meanFnBin = meanBin/nCMULBin;
3908 Double_t meanFnBinError = TMath::Sqrt( TMath::Power(meanBinError/nCMULBin,2) + TMath::Power(meanBin/TMath::Power(nCMULBin,2.),2) );
3909
3910 if (printout) std::cout << meanBinSys/nCMULBin/meanFnBin << std::endl;
a58729a5 3911
cb690a12 3912 if (printout) std::cout << meanFnBin << " +- " << meanFnBinError << std::endl;
a58729a5 3913
cb690a12 3914 hFnormMean->SetBinContent(i,meanFnBin);
3915 hFnormMean->SetBinError(i,meanFnBinError);
a58729a5 3916 }
2754fd07 3917
cb690a12 3918 TH1* o = fMergeableCollection->Histo(id.Data(),hMBMean->GetName());
a58729a5 3919
cb690a12 3920 if (o)
2754fd07 3921 {
cb690a12 3922 AliWarning(Form("Replacing %s/%s",id.Data(),hMBMean->GetName()));
3923 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hMBMean->GetName()));
2754fd07 3924 }
2f331ac9 3925
cb690a12 3926 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hMBMean);
a58729a5 3927
cb690a12 3928 if ( adoptOK ) std::cout << "+++FNorm histo " << hMBMean->GetName() << " adopted" << std::endl;
3929 else AliError(Form("Could not adopt FNorm histo %s",hMBMean->GetName()));
a58729a5 3930
a58729a5 3931
cb690a12 3932 o = fMergeableCollection->Histo(id.Data(),hFnormMean->GetName());
a58729a5 3933
cb690a12 3934 if (o)
a58729a5 3935 {
cb690a12 3936 AliWarning(Form("Replacing %s/%s",id.Data(),hFnormMean->GetName()));
3937 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFnormMean->GetName()));
a58729a5 3938 }
2f331ac9 3939
cb690a12 3940 adoptOK = fMergeableCollection->Adopt(id.Data(),hFnormMean);
a58729a5 3941
cb690a12 3942 if ( adoptOK ) std::cout << "+++FNorm histo " << hFnormMean->GetName() << " adopted" << std::endl;
3943 else AliError(Form("Could not adopt FNorm histo %s",hFnormMean->GetName()));
3944
3945
2f331ac9 3946}
3947
3948//_____________________________________________________________________________
cb690a12 3949void AliAnalysisMuMu::ComputeIntFnormFromCounters(const char* triggerCluster, const char* eventSelection, const char* filePileUpCorr, Bool_t printout)
2f331ac9 3950{
cb690a12 3951 /// Compute the CMUL to CINT ratio(s) in 2 steps from the CC(), integrated
2f331ac9 3952
cb690a12 3953 TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3954 if ( colType.Contains("-B-") ) colType = "B";
3955 else if ( colType.Contains("-S-") ) colType = "S";
3956 else
2f331ac9 3957 {
cb690a12 3958 AliError("Unknown collision type");
3959 return;
2f331ac9 3960 }
2f331ac9 3961
cb690a12 3962 TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3963 if ( triggerType.Contains("7-") ) triggerType = "7";
3964 else if ( triggerType.Contains("8-") ) triggerType = "8";
3965 else
a58729a5 3966 {
cb690a12 3967 AliError("Unknown trigger type");
3968 return;
a58729a5 3969 }
cb690a12 3970
3971 TString striggerCluster(triggerCluster);
3972 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3973 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3974 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3975 else
a58729a5 3976 {
cb690a12 3977 AliError("Unknown trigger cluster");
3978 return;
a58729a5 3979 }
cb690a12 3980
3981 TString seventSelection(eventSelection);
3982 TString sruns = CC()->GetKeyWords("run");
3983 TObjArray* runs = sruns.Tokenize(",");
3984 Double_t NofRuns = runs->GetEntries();
2f331ac9 3985
cb690a12 3986 Bool_t corrPU(kFALSE);
3987 TObjArray* pUCorr = new TObjArray();
3988 if ( strlen(filePileUpCorr) > 0 )
a58729a5 3989 {
cb690a12 3990// std::cout << "Extracting Pile-Up correction factors from " << filePileUpCorr << std::endl;
3991 char line[1024];
3992 ifstream in(filePileUpCorr);
a58729a5 3993
cb690a12 3994 while ( in.getline(line,1024,'\n'))
a58729a5 3995 {
cb690a12 3996 TString lrun(line);
3997 TString lvalue(line);
3998
3999 lrun.Remove(0,4);
4000 lrun.Remove(6,67);
4001
4002 lvalue.Remove(0,57);//71
4003
4004// std::cout << "RUN: " << lrun.Data() << " PUFactor = " << lvalue.Data() << std::endl;
4005
4006 pUCorr->Add(new TParameter<Double_t>(lrun.Data(),lvalue.Atof()));
a58729a5 4007 }
cb690a12 4008 corrPU = kTRUE;
a58729a5 4009 }
a58729a5 4010
cb690a12 4011 TIter nextRun(runs);
4012 TObjString* s;
4013
4014 TH1* h;
a58729a5 4015
cb690a12 4016 Double_t FNormTot(0.);
4017 Double_t FNormTotError(0.);
2f331ac9 4018
cb690a12 4019 TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
2f331ac9 4020
cb690a12 4021 h = new TH1F("hFNormIntVSrun","Integrated Normalization factor vs run;run;FNorm",NofRuns,1.,NofRuns);
2f331ac9 4022
cb690a12 4023 Double_t nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
4024 seventSelection.Data(),triggerType.Data(),colType.Data())); //Total Nof CMUL7 events
2f331ac9 4025
cb690a12 4026 TObjArray* aCluster = striggerCluster.Tokenize("-");
4027 TIter nextCluster(aCluster);
4028 TObjString* striggerClusterS;
2f331ac9 4029
cb690a12 4030 TList* lRun2Reject = new TList();
4031 lRun2Reject->SetOwner(kTRUE);
2f331ac9 4032
cb690a12 4033 nextRun.Reset();
4034 Int_t i(0);//Run label index
4035 std::cout << std::endl;
4036 std::cout << std::endl;
4037 std::cout << std::endl;
4038 std::cout << "Computing FNorm" << std::endl;
4039 while ( ( s = static_cast<TObjString*>(nextRun())) ) //Run loop
a58729a5 4040 {
cb690a12 4041 Double_t nCMSL(0.),nCMSLandOMUL(0.);
4042 nextCluster.Reset();
4043 while ( (striggerClusterS = static_cast<TObjString*>(nextCluster())) && nCMSL == 0. )
4044 {
4045 nCMSL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s/centrality:V0A/run:%s",
4046 seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName()));
4047
4048 nCMSLandOMUL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s&0MUL/centrality:V0A/run:%s",
4049 seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName()));
4050 }
81190958 4051
cb690a12 4052 Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/run:%s",
4053 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
a58729a5 4054
cb690a12 4055 Double_t nCINTandOMSL = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD&0MSL/centrality:V0A/run:%s",
4056 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
a58729a5 4057
cb690a12 4058 Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s",
4059 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
81190958 4060
cb690a12 4061 Double_t FNorm(0.),FNormError(0.);
4062 Double_t pUfactor = 1.;
4063 if ( nCMSLandOMUL != 0. && nCINTandOMSL !=0. && nCMSL != 0. && nCINT !=0. )
4064 {
4065 if (corrPU)
4066 {
4067 TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(pUCorr->FindObject(s->GetName()));
4068 if ( p ) pUfactor = p->GetVal();
4069 else
4070 {
4071 AliError(Form("Run %s not found in pile-up correction list",s->GetName()));
4072 }
4073 }
4074 FNorm = (nCMSL*nCINT)*pUfactor/(nCMSLandOMUL*nCINTandOMSL);
4075 FNormError = ErrorPropagationAxBoverCxD(nCMSL,nCINT,nCMSLandOMUL,nCINTandOMSL)*pUfactor;
4076 }
4077 else
81190958 4078 {
cb690a12 4079 if ( nCINT == 0 ) std::cout << " Warning: Bad run " << s->GetName() << " has no MB trigger in this bin. Remove from analysis" << std::endl;
4080
4081 lRun2Reject->Add(new TObjString(s->GetName()));
4082 if ( printout ) std::cout << "Run " << s->GetName() << " not used for FNorm cause lack of stats" << std::endl;
4083 continue;
81190958 4084 }
cb690a12 4085
4086 FNormTot += FNorm*nCMUL; // This is the sum of equivalent Nof MB per CMUL run by run. NOTE: This sum is NOT always the total equivalent Nof MB per CMUL because in pp 2012 if just one cluster is used at a time this sum is not the sum for all runs
4087 FNormTotError += TMath::Power(nCMUL*FNormError,2.) + TMath::Power(FNorm*TMath::Sqrt(nCMUL),2.);
4088
4089 if ( printout ) std::cout << "Run " << s->GetName() << " FNorm = " << FNorm << " +- " << FNormError << " ; PUFactor =" << pUfactor << " ; " << "Nof CMUL = " << nCMUL << std::endl;
4090
4091 h->GetXaxis()->SetBinLabel(++i,s->GetName());
4092 h->SetBinContent(i,FNorm);
4093 h->SetBinError(i,FNormError);
4094
a58729a5 4095 }
2f331ac9 4096
cb690a12 4097 TIter nextRejectRun(lRun2Reject);
4098 TObjString* run2Rej;
4099 Double_t nCMULTotRej(0.);
4100 while ( (run2Rej = static_cast<TObjString*>(nextRejectRun())) )
4101 {
4102 nCMULTotRej += CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s",
4103 seventSelection.Data(),triggerType.Data(),colType.Data(),
4104 run2Rej->GetName())); //Sum of CMUL7 events from rejected runs
4105 }
2f331ac9 4106
cb690a12 4107 nCMULTot = nCMULTot - nCMULTotRej;
4108
4109 TH1* o = fMergeableCollection->Histo(id.Data(),h->GetName());
4110
4111 if (o)
2f331ac9 4112 {
cb690a12 4113 AliWarning(Form("Replacing %s/%s",id.Data(),h->GetName()));
4114 fMergeableCollection->Remove(Form("%s/%s",id.Data(),h->GetName()));
2f331ac9 4115 }
2f331ac9 4116
cb690a12 4117 Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),h);
a58729a5 4118
cb690a12 4119 if ( adoptOK ) std::cout << "+++FNorm histo " << h->GetName() << " adopted" << std::endl;
4120 else AliError(Form("Could not adopt FNorm histo %s",h->GetName()));
4121
4122 //___
4123
4124 FNormTotError = TMath::Sqrt(TMath::Power(TMath::Sqrt(FNormTotError)/nCMULTot,2.) + TMath::Power(FNormTot*TMath::Sqrt(nCMULTot)/TMath::Power(nCMULTot,2.),2.));
4125
4126 FNormTot = FNormTot/nCMULTot; // nCMULTot is here nCMULTot - nCMULTotRej
4127
4128 std::cout << "FNorm = " << FNormTot << " +- " << FNormTotError << std::endl;
4129
4130 TH1* hFNormTot = new TH1F("hFNormInt","Global Normalization factor",1,0.,1.);
4131 hFNormTot->SetBinContent(1,FNormTot);
4132 hFNormTot->SetBinError(1,FNormTotError);
4133
4134 o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
4135
4136 if (o)
2f331ac9 4137 {
cb690a12 4138 AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
4139 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
2f331ac9 4140 }
4141
cb690a12 4142 adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
a58729a5 4143
cb690a12 4144 if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
4145 else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
4146
4147 //___
4148 TH1* hNEqMB = new TH1F("hNEqMB","Equivalent number of MB events per CMUL",1,0.,1.);
4149
4150 nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
4151 seventSelection.Data(),triggerType.Data(),colType.Data())); //Total Nof CMUL7 events
4152
4153 Double_t nofEqMB = FNormTot*nCMULTot;
4154 Double_t nofEqMBError = TMath::Sqrt( TMath::Power(FNormTotError*nCMULTot,2.) + TMath::Power(FNormTot*TMath::Sqrt(nCMULTot),2.) );
4155
4156 std::cout << "EqMB = " << nofEqMB << " +- " << TMath::Sqrt(nofEqMBError) << std::endl;
4157
4158 hNEqMB->SetBinContent(1,nofEqMB);
4159 hNEqMB->SetBinError(1,nofEqMBError);
4160
4161 o = fMergeableCollection->Histo(id.Data(),hNEqMB->GetName());
4162
4163 if (o)
a58729a5 4164 {
cb690a12 4165 AliWarning(Form("Replacing %s/%s",id.Data(),hNEqMB->GetName()));
4166 fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNEqMB->GetName()));
a58729a5 4167 }
4168
cb690a12 4169 adoptOK = fMergeableCollection->Adopt(id.Data(),hNEqMB);
2f331ac9 4170
cb690a12 4171 if ( adoptOK ) std::cout << "+++FNorm histo " << hNEqMB->GetName() << " adopted" << std::endl;
4172 else AliError(Form("Could not adopt FNorm histo %s",hNEqMB->GetName()));
2f331ac9 4173
cb690a12 4174 //___
2f331ac9 4175
cb690a12 4176 delete runs;
4177 delete lRun2Reject;
4178 delete aCluster;
4179
4180 return;
2f331ac9 4181
2f331ac9 4182}
4183
4184//_____________________________________________________________________________
cb690a12 4185void AliAnalysisMuMu::PlotYiedWSyst(const char* triggerCluster)
4186{
4187 TString striggerCluster(triggerCluster);
4188 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
4189 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
4190 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
4191 else
2f331ac9 4192 {
cb690a12 4193 AliError("Unknown trigger cluster");
4194 return;
2f331ac9 4195 }
4196
cb690a12 4197 TString path(Form("%s/%s/%s/%s",
4198 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4199 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4200 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4201 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
4202
4203 TH1* hY = OC()->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),"hJPsiYieldVSdNchdEtaRelative");
4204 if ( !hY )
1afce1ce 4205 {
cb690a12 4206 AliError("No yield found");
4207 return;
1afce1ce 4208 }
4209
cb690a12 4210 TString id(Form("/TESTSYST/%s",path.Data()));
1afce1ce 4211
cb690a12 4212 TH1* hYSyst = static_cast<TH1*>(hY->Clone("RelYieldSyst"));
4213 if ( !hYSyst )
1afce1ce 4214 {
cb690a12 4215 AliError("No systematic found");
4216 return;
1afce1ce 4217 }
a58729a5 4218
cb690a12 4219 TH1* hS = OC()->Histo(id.Data(),"yield_Systematics");
1afce1ce 4220
cb690a12 4221 for ( Int_t i = 1 ; i <= hY->GetNbinsX() ; i++ )
4222 {
4223 hYSyst->SetBinError(i,hS->GetBinContent(i)*hY->GetBinContent(i)/100.);
4224 }
4225
4226 hY->Draw();
4227 hYSyst->Draw("same");
a58729a5 4228}
4229
cb690a12 4230
a58729a5 4231//_____________________________________________________________________________
cb690a12 4232void AliAnalysisMuMu::ComputeJpsiYield(AliMergeableCollection* oc, Bool_t relative, const char* fNormType, const char* triggerCluster,
4233 const char* whatever, const char* sResName, AliMergeableCollection* ocMBTrigger, Double_t mNTrCorrection)
a58729a5 4234{
cb690a12 4235 // This method is suppossed to be used from the file with the counters, oc is the AliMergeableCollection of the file with the histograms (if separated, which is better since we do not need the minv,mean pt... analysis in CINT&0MUL... triggers)
4236 // ocMBTrigger is the mergeableCollection with the MB trigger dNchdEta plot (migth be the same as oc, in which case we set ocMBTrigger=0x0)
4237 //FIXME::Make it general
4238
4239 TString sfNormType(fNormType);
4240 TString swhat("");
4241 TString sres("");
4242 TString swhatever(whatever);
4243 if ( swhatever.Contains("DNCHDETA"))
4244 {
4245 swhat = "dNchdEta";
4246 if ( strlen(sResName) > 0 ) sres = sResName; //"PSIPSIPRIMECB2VWGINDEPTAILS";
4247 }
4248 else if ( swhatever.Contains("NTRCORR") )
4249 {
4250 swhat = "Nch";
4251 if ( strlen(sResName) > 0 ) sres = sResName; //"PSIPSIPRIMECB2VWG_2.0_5.0";
4252 }
2f331ac9 4253
cb690a12 4254 if ( IsSimulation() )
a58729a5 4255 {
cb690a12 4256 AliError("Cannot compute J/Psi yield: Is a simulation file");
4257 return;
a58729a5 4258 }
cb690a12 4259
4260 TString striggerCluster(triggerCluster);
4261 if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
4262 else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
4263 else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
a58729a5 4264 else
4265 {
cb690a12 4266 AliError("Unknown trigger cluster");
4267 return;
a58729a5 4268 }
1afce1ce 4269
cb690a12 4270 TString path(Form("%s/%s/%s/%s",
4271 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4272 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4273 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4274 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
2f331ac9 4275
cb690a12 4276 AliMergeableCollection* mc;
4277 if ( !oc ) mc = OC();
4278 else mc = oc;
2f331ac9 4279
cb690a12 4280 Double_t bR = 0.0593; // BR(JPsi->mu+mu-)
4281 Double_t bRerror = 0.0006 ;
2f331ac9 4282
cb690a12 4283 //_________Integrated yield
4284 AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(mc->GetObject(Form("/%s/%s",path.Data(),"PSI-INTEGRATED-AccEffCorr")));
4285 if ( !sInt )
4286 {
4287 AliError(Form("No spectra %s found in %s","PSI-INTEGRATED-AccEffCorr",path.Data()));
4288 return;
4289 }
2f331ac9 4290
cb690a12 4291 AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
4292 b->AddBin("psi","INTEGRATED");
2f331ac9 4293
cb690a12 4294 AliAnalysisMuMuBinning::Range* bin = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
2f331ac9 4295
cb690a12 4296 AliAnalysisMuMuResult* result = sInt->GetResultForBin(*bin);
4297 if ( !result )
4298 {
4299 AliError(Form("No result for bin %s found in %s",bin->AsString().Data(),"PSI-INTEGRATED-AccEffCorr"));
4300 return;
4301 }
2f331ac9 4302
cb690a12 4303// if ( strlen(sResName) > 0/*sResName.Sizeof() > 0*/ )
4304// {
4305// result = result->SubResult(sres.Data());//INDEPTAILS
4306// if ( !result )
4307// {
4308// AliError(Form("No subresult %s found in %s",sres.Data(),path.Data()));
4309// return;
4310// }
4311// }
2f331ac9 4312
cb690a12 4313 Double_t NofJPsiTot = result->GetValue("NofJPsi");
4314 Double_t NofJPsiTotError = result->GetErrorStat("NofJPsi");
2f331ac9 4315
cb690a12 4316 TH1* hMBTot = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNEqMB",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4317 if ( !hMBTot )
4318 {
4319 AliError(Form("No eq Nof MB events found in %s",Form("/FNORM-%s/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/V0A/hNEqMB",striggerCluster.Data())));
4320 return;
4321 }
2f331ac9 4322
cb690a12 4323 Double_t nEqMBTot = hMBTot->GetBinContent(1);
4324 Double_t nEqMBTotError = hMBTot->GetBinError(1);
4325
4326 Double_t yieldInt = NofJPsiTot/(nEqMBTot*bR);
4327 Double_t yieldIntError = TMath::Sqrt(TMath::Power(NofJPsiTotError/(nEqMBTot*bR),2.) +
4328 TMath::Power(nEqMBTotError*NofJPsiTot*bR/TMath::Power(nEqMBTot*bR,2.),2.) +
4329 TMath::Power(NofJPsiTot*nEqMBTot*bRerror/TMath::Power(nEqMBTot*bR,2.),2.));
4330
4331 std::cout << "Integrated yield = " << yieldInt << " +- " << yieldIntError << std::endl;
4332
4333 TH1* hYint = new TH1F("hJPsiYieldInt","Integrated J/#psi yield",1,0.,1.);
4334 hYint->SetBinContent(1,yieldInt);
4335 hYint->SetBinError(1,yieldIntError);
4336
4337 TH1* o = mc->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hYint->GetName());
4338
4339 if (o)
2f331ac9 4340 {
cb690a12 4341 AliWarning(Form("Replacing /RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hYint->GetName()));
4342 mc->Remove(Form("/RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hYint->GetName()));
4343 }
4344
4345 Bool_t adoptOK = mc->Adopt(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hYint);
4346
4347 if ( adoptOK ) std::cout << "+++Yield histo " << hYint->GetName() << " adopted" << std::endl;
4348 else AliError(Form("Could not adopt Yield histo %s",hYint->GetName()));
4349
4350 delete b;
4351
4352 //_____Differential yield
4353
4354 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(mc->GetObject(Form("/%s/%s",path.Data(),whatever)));
4355 if ( !s )
4356 {
4357 AliError(Form("No spectra %s found in %s",whatever,path.Data()));
4358 return;
4359 }
4360
4361 std::cout << "Number of J/Psi:" << std::endl;
4362 TH1* hry = s->Plot("NofJPsi",sres.Data(),kFALSE);//INDEPTAILS //Number of Jpsi
4363
4364 std::cout << "" << std::endl;
4365
4366// std::cout << "Equivalent number of MB events:" << std::endl;
4367 TH1* hMB(0x0);
4368 if ( sfNormType.Contains("offline") )
4369 {
4370 hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEta",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4371 if ( !hMB )
a58729a5 4372 {
cb690a12 4373 AliError("Histo hNofEqMBVSdNchdEta not found");
4374 return;
a58729a5 4375 }
cb690a12 4376 }
4377 else if ( sfNormType.Contains("global") )
4378 {
4379 hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEtaFromGlobal",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4380 if ( !hMB )
2f331ac9 4381 {
cb690a12 4382 AliError("Histo hNofEqMBVSdNchdEtaFromGlobal not found");
4383 return;
a58729a5 4384 }
cb690a12 4385 }
4386 else if ( sfNormType.Contains("mean") )
4387 {
4388 hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEtaFromMean",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4389 if ( !hMB )
a58729a5 4390 {
cb690a12 4391 AliError("Histo hNofEqMBVSdNchdEtaFromMean not found");
4392 return;
2f331ac9 4393 }
2f331ac9 4394 }
cb690a12 4395 else
4396 {
4397 AliError("Dont know what Fnorm use");
4398 return;
4399 }
2f331ac9 4400
cb690a12 4401 TH1* hy;
4402 if ( relative )
a58729a5 4403 {
cb690a12 4404 TString path2(Form("/%s/%s/%s",
4405 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4406 First(Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,kFALSE)).Data(),
4407 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
4408
4409 TH1* hdNch;
4410 if ( ocMBTrigger ) hdNch = ocMBTrigger->Histo(path2.Data(),swhat.Data());//dNchdEta
4411 else hdNch = mc->Histo(path2.Data(),swhat.Data());//dNchdEta
4412
4413 const TArrayD* binArray = hry->GetXaxis()->GetXbins();
4414 Int_t size = binArray->GetSize();
4415 Double_t* axis = new Double_t[size];
4416 for ( Int_t k = 0 ; k < size ; k++ )
0de0b867 4417 {
cb690a12 4418 axis[k] = binArray->At(k)/(hdNch->GetMean()*(1 - mNTrCorrection));
0de0b867 4419 }
cb690a12 4420
4421 hy = new TH1D("hJPsiYieldVSdNchdEtaRelative","Relative J/#psi yield vs dN_{ch}/d#eta;dN_{ch}/d#eta/<dN_{ch}/d#eta>;Y^{J/#psi}/Y^{J/#psi}_{int}",size-1,axis);
4422 delete axis;
a58729a5 4423 }
cb690a12 4424 else
4425 {
4426 hy = static_cast<TH1D*>(hry->Clone("hJPsiYieldVSdNchdEta"));
4427 hy->SetTitle("J/#psi yield vs dN_{ch}/d#eta");
4428 hy->GetXaxis()->SetTitle("dN_{ch}/d#eta");
4429 hy->GetYaxis()->SetTitle("Y^{J/#psi}");
4430 }
4431 // AccxEff(from rel diff or paper) // Signal extraction
4432// Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.01,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.022,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.009,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.016 ,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.033,2.) )}; //FIXME: find a way to give this as input
4433// Double_t systFNorm[9] = {0.003,0.001,0.002,0.003,0.002,0.004,0.011,0.012,0.071};
4434// Double_t systPU[9] = {0.00,0.01,0.012,0.014,0.014,0.019,0.020,0.021,0.040}; //_______pPb
4435 // AccxEff(from paper)
4436// Double_t systNofJpsiTot = 0.015;
4437
4438// Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.028,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.016,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.024,2.) + TMath::Power(0.005 ,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.016,2.) )}; //FIXME: find a way to give this as input
4439// Double_t systFNorm[9] = {0.005,0.004,0.004,0.004,0.003,0.002,0.002,0.04,0.04};
4440// Double_t systPU[9] = {0.00,0.007,0.015,0.011,0.014,0.018,0.014,0.011,0.020}; //______Pbp
4441// Double_t systNofJpsiTot = 0.015;
4442
4443// Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.034,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.042,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.063,2.) + TMath::Power(0.014,2.) ),TMath::Sqrt( TMath::Power(0.094,2.) + TMath::Power(0.009,2.) ),TMath::Sqrt( TMath::Power(0.00,2.) + TMath::Power(0.00 ,2.) ),TMath::Sqrt( TMath::Power(0.00,2.) + TMath::Power(0.00,2.) )}; //FIXME: find a way to give this as input
4444// Double_t systFNorm[9] = {0.004,0.019,0.002,0.012,0.048,0.063,0.082,0.000,0.000};
4445// Double_t systPU[9] = {0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00}; //______pp |eta|<0.5
4446// Double_t systNofJpsiTot = 0.017;
4447
4448 Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.037,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.021,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.022,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.019,2.) + TMath::Power(0.001,2.) ),TMath::Sqrt( TMath::Power(0.036,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.042,2.) + TMath::Power(0.001,2.) ),TMath::Sqrt( TMath::Power(0.039,2.) + TMath::Power(0.012 ,2.) ),TMath::Sqrt( TMath::Power(0.000,2.) + TMath::Power(0.000,2.) )}; //FIXME: find a way to give this as input
4449 Double_t systFNorm[9] = {0.026,0.002,0.015,0.019,0.012,0.030,0.015,0.119,0.000};
4450 Double_t systPU[9] = {0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00}; //______pp |eta|<1
4451 Double_t systNofJpsiTot = 0.017;
4452
4453 for ( Int_t i = 1 ; i <= hy->GetNbinsX() ; i++ )
4454 {
4455 Double_t yield = hry->GetBinContent(i)/(hMB->GetBinContent(i)*bR);
4456 Double_t yieldError = TMath::Sqrt(TMath::Power(hry->GetBinError(i)/(hMB->GetBinContent(i)*bR),2.) +
4457 TMath::Power(hMB->GetBinError(i)*hry->GetBinContent(i)*bR/TMath::Power(hMB->GetBinContent(i)*bR,2.),2.) +
4458 TMath::Power(hry->GetBinContent(i)*hMB->GetBinContent(i)*bRerror/TMath::Power(hMB->GetBinContent(i)*bR,2.),2.));
4459
4460// std::cout << "Differential yield bin " << i << " = " << yield << " +- " << yieldError << std::endl;
4461
4462 if ( relative )
4463 {
4464 yieldError = TMath::Sqrt(TMath::Power(yieldError/yieldInt,2.) + TMath::Power((yield*yieldIntError)/TMath::Power(yieldInt,2.),2.));
4465 yield /= yieldInt;
4466
4467// std::cout << "relative yield bin " << i << " = " << yield << " +- " << yieldError << std::endl;
4468 Double_t sNJpsiBin = hry->GetBinContent(i)*systNofJpsiBin[i-1];
4469 Double_t sNJpsiTot = NofJPsiTot*systNofJpsiTot;
4470 Double_t sMBBin = hMB->GetBinContent(i)*systFNorm[i-1];
4471 Double_t sMBTot = nEqMBTot*0.01;
4472
4473 Double_t syst = TMath::Sqrt( TMath::Power((sNJpsiBin/NofJPsiTot)*(nEqMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)*sNJpsiTot/TMath::Power(NofJPsiTot,2.))*(nEqMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)/NofJPsiTot)*(sMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)/NofJPsiTot)*(sMBBin*nEqMBTot/TMath::Power(hMB->GetBinContent(i),2.)),2.) );
4474
4475 std::cout << "sys" << syst/yield << " w/pu = " << TMath::Sqrt( TMath::Power(syst/yield,2.) + TMath::Power(systPU[i-1],2.)) << std::endl;
4476 std::cout << yield << " +- " << yieldError << std::endl;
4477 }
a58729a5 4478
cb690a12 4479 hy->SetBinContent(i,yield);
4480 hy->SetBinError(i,yieldError);
4481 }
a58729a5 4482
cb690a12 4483 o = mc->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hy->GetName());
81190958 4484
cb690a12 4485 if (o)
81190958 4486 {
cb690a12 4487 AliWarning(Form("Replacing %s/%s","/RESULTS-%s/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/V0A",hy->GetName()));
4488 mc->Remove(Form("/RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hy->GetName()));
81190958 4489 }
cb690a12 4490
4491 adoptOK = mc->Adopt(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hy);
4492
4493 if ( adoptOK ) std::cout << "+++Yield histo " << hy->GetName() << " adopted" << std::endl;
4494 else AliError(Form("Could not adopt Yield histo %s",hy->GetName()));
1afce1ce 4495
81190958 4496
cb690a12 4497
4498 delete hry;
1afce1ce 4499
a58729a5 4500
cb690a12 4501 return;
a58729a5 4502}
4503
4504//_____________________________________________________________________________
cb690a12 4505void AliAnalysisMuMu::ComputeJpsiMPt(Bool_t relative, const char* whatever, const char* sResName, AliMergeableCollection* ocMBTrigger, Double_t mNTrCorrection)
a58729a5 4506{
cb690a12 4507 // ocMBTrigger is the mergeableCollection with the MB trigger dNchdEta plot (migth be the same as oc, in which case we set ocMBTrigger=0x0)
4508 //FIXME::Make it general
4509
4510
4511 TString swhat("");
4512 TString sres("");
4513 TString swhatever(whatever);
4514// if ( swhatever.Contains("DNCHDETA"))
4515// {
4516// swhat = "dNchdEta";
4517// sres = "MPT2CB2VWGPOL2INDEPTAILS";
4518// }
4519// else
4520 if ( swhatever.Contains("NTRCORR") )
4521 {
4522 swhat = "Nch";
4523 if ( strlen(sResName) > 0 ) sres = sResName; //sres = "MPTPSIPSIPRIMECB2VWG_BKGMPTPOL2";
4524 }
a58729a5 4525
cb690a12 4526 if ( IsSimulation() )
4527 {
4528 AliError("Cannot compute J/Psi yield: Is a simulation file");
4529 return;
4530 }
4531
4532 TString path(Form("%s/%s/%s/%s",
4533 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4534 First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4535 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4536 First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
4537
4538 //_________Integrated mean pt
4539 AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),"PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS")));
4540 if ( !sInt )
4541 {
4542 AliError(Form("No spectra %s found in %s","PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS",path.Data()));
4543 return;
4544 }
4545
4546 AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
4547 b->AddBin("psi","INTEGRATED");
4548
4549 AliAnalysisMuMuBinning::Range* bin = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
4550
4551 AliAnalysisMuMuResult* result = sInt->GetResultForBin(*bin);
4552 if ( !result )
4553 {
4554 AliError(Form("No result for bin %s found in spectra %s",bin->AsString().Data(),sInt->GetName()));
4555 return;
4556 }
4557
4558// if ( sres.Sizeof() > 0 )
4559// {
4560// result = result->SubResult(sres.Data());
4561// if ( !result )
4562// {
4563// AliError(Form("No subresult %s found in result",result->GetName()));
4564// return;
4565// }
4566// AliAnalysisMuMuResult* subresult = result->SubResult(sres.Data());//"MPT2CB2VWGPOL2INDEPTAILS"
4567// if ( !subresult )
4568// {
4569// AliError(Form("No subresult MPT2CB2VWGPOL2 found in result %s",result->GetName()));
4570// return;
4571// }
4572
4573// }
4574// Double_t JPsiMPtTot = subresult->GetValue("MeanPtJPsi");
4575// Double_t JPsiMPtTotError = subresult->GetErrorStat("MeanPtJPsi");
4576
4577 Double_t JPsiMPtTot = result->GetValue("MeanPtJPsi");
4578 Double_t JPsiMPtTotError = result->GetErrorStat("MeanPtJPsi");
4579
4580 TH1* hMPtint = new TH1F("hJPsiMPtInt","Integrated J/#psi mean p_{T}",1,0.,1.);
4581 hMPtint->SetBinContent(1,JPsiMPtTot);
4582 hMPtint->SetBinError(1,JPsiMPtTotError);
a58729a5 4583
cb690a12 4584 TH1* o = OC()->Histo(Form("/RESULTS/%s",path.Data()),hMPtint->GetName());
a58729a5 4585
cb690a12 4586 if (o)
a58729a5 4587 {
cb690a12 4588 AliWarning(Form("Replacing /RESULTS/%s/%s",path.Data(),hMPtint->GetName()));
4589 OC()->Remove(Form("/RESULTS/%s/%s",path.Data(),hMPtint->GetName()));
4590 }
a58729a5 4591
cb690a12 4592 Bool_t adoptOK = OC()->Adopt(Form("/RESULTS/%s",path.Data()),hMPtint);
4593
4594 if ( adoptOK ) std::cout << "+++Mean Pt histo " << hMPtint->GetName() << " adopted" << std::endl;
4595 else AliError(Form("Could not adopt Mean Pt histo %s",hMPtint->GetName()));
4596
4597 delete b;
4598
4599 //_____Differential mean pt
4600
4601 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),whatever)));
4602 if ( !s )
4603 {
4604 AliError(Form("No spectra %s found in %s",whatever,path.Data()));
4605 return;
4606 }
4607
4608 std::cout << "Mean pt of J/Psi:" << std::endl;
4609 TH1* hrmPt = s->Plot("MeanPtJPsi",sres.Data(),kFALSE); //MPT2CB2VWGPOL2INDEPTAILS//mean pt of Jpsi
4610 std::cout << "" << std::endl;
4611
4612 Double_t ptInt,ptIntError;
4613 TH1* hmPt;
4614 if ( relative )
4615 {
4616 TString path2(Form("/%s/%s/%s",
4617 First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4618 First(Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,kFALSE)).Data(),
4619 First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
4620
4621 TH1* hdNch;
4622 if ( ocMBTrigger ) hdNch = ocMBTrigger->Histo(path2.Data(),swhat.Data());
4623 else hdNch = OC()->Histo(path2.Data(),swhat.Data());
4624
4625 const TArrayD* binArray = hrmPt->GetXaxis()->GetXbins();
4626 Int_t size = binArray->GetSize();
4627 Double_t* axis = new Double_t[size];
4628 for ( Int_t k = 0 ; k < size ; k++ )
4629 {
4630 axis[k] = binArray->At(k)/(hdNch->GetMean()*(1 - mNTrCorrection));
4631 }
a58729a5 4632
cb690a12 4633 hmPt = new TH1D("hJPsiMeanPtVSdNchdEtaRelative","Relative J/#psi mean p_{T} vs dN_{ch}/d#eta/<dN_{ch}/d#eta>;dN_{ch}/d#eta/<dN_{ch}/d#eta>;<p_{T}^{J/#psi}>/<p_{T}^{J/#psi}_{int}>",size-1,axis);
4634 delete axis;
a58729a5 4635
cb690a12 4636 ptInt = result->GetValue("MeanPtJPsi",sres.Data());
4637 ptIntError = result->GetErrorStat("MeanPtJPsi",sres.Data());
4638
4639// delete b;
a58729a5 4640 }
4641 else
4642 {
cb690a12 4643 hmPt = static_cast<TH1D*>(hrmPt->Clone("hJPsiMeanPtVSdNchdEta"));
4644 hmPt->SetTitle("J/#psi mean p_{T} vs dN_{ch}/d#eta");
4645 hmPt->GetXaxis()->SetTitle("dN_{ch}/d#eta");
4646 hmPt->GetYaxis()->SetTitle("<p_{T}^{J/#psi}>");
4647 }
4648
4649 Double_t systMptInt[9] = {0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014}; //FIXME: find a way to give this as input
4650 Double_t systMptBin[9] = {0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014};
4651
4652 Double_t systMptRel[9] = {0.002,0.001,0.001,0.002,0.002,0.002,0.002,0.004,0.004}; //signal extraction pPb
4653
4654// Double_t systMptRel[9] = {0.002,0.001,0.012,0.001,0.002,0.002,0.001,0.004,0.003}; //signal extraction Pbp
4655
4656// Double_t systMptRel[9] = {0.001,0.002,0.001,0.002,0.002,0.003,0.005,0.000,0.000}; //signal extraction pp|eta|<05
4657
4658// Double_t systMptRel[9] = {0.002,0.002,0.002,0.002,0.001,0.002,0.001,0.012,0.000}; //signal extraction pp|eta|<1
4659
4660 for ( Int_t i = 1 ; i <= hrmPt->GetNbinsX() ; i++ )
4661 {
4662 Double_t pt = hrmPt->GetBinContent(i);
4663 Double_t ptError = hrmPt->GetBinError(i);
4664
4665 if ( relative )
a58729a5 4666 {
cb690a12 4667 ptError = TMath::Sqrt(TMath::Power(ptError/ptInt,2.) + TMath::Power((pt*ptIntError)/TMath::Power(ptInt,2.),2.));
a58729a5 4668
cb690a12 4669 Double_t sMptInt = ptInt*systMptInt[i-1];
4670 Double_t sMptBin = pt*systMptBin[i-1];
4671 Double_t sysMptRel = TMath::Sqrt( TMath::Power(sMptBin/ptInt,2) + TMath::Power(pt*sMptInt/TMath::Power(ptInt,2.),2.) );
a58729a5 4672
cb690a12 4673 pt /= ptInt;
4674
4675 std::cout << TMath::Sqrt( TMath::Power(sysMptRel/pt,2.) +TMath::Power(systMptRel[i-1],2.) ) << std::endl;
4676
4677 std::cout << pt << " +- " << ptError << std::endl;
4678
a58729a5 4679 }
cb690a12 4680
4681 hmPt->SetBinContent(i,pt);
4682 hmPt->SetBinError(i,ptError);
4683 }
4684
4685 o = fMergeableCollection->Histo(Form("/RESULTS/%s",path.Data()),hmPt->GetName());
4686
4687 if (o)
4688 {
4689 AliWarning(Form("Replacing /RESULTS/%s/%s",path.Data(),hmPt->GetName()));
4690 fMergeableCollection->Remove(Form("/RESULTS/%s/%s",path.Data(),hmPt->GetName()));
a58729a5 4691 }
cb690a12 4692
4693 adoptOK = fMergeableCollection->Adopt(Form("/RESULTS/%s",path.Data()),hmPt);
4694
4695 if ( adoptOK ) std::cout << "+++Mean Pt histo " << hmPt->GetName() << " adopted" << std::endl;
4696 else AliError(Form("Could not adopt mean pt histo %s",hmPt->GetName()));
4697
4698
4699
4700 delete hrmPt;
4701
4702
4703 return;
a58729a5 4704
a58729a5 4705
a58729a5 4706}
4707
81190958 4708//_____________________________________________________________________________
cb690a12 4709Double_t AliAnalysisMuMu::ErrorPropagationAxBoverCxD(Double_t a,Double_t b,Double_t c,Double_t d)
81190958 4710{
cb690a12 4711 //Just valid for counts
4712 Double_t error2 = TMath::Power(b/(c*d),2.)*a + TMath::Power(a/(c*d),2.)*b + TMath::Power(a*b*d,2.)*(c/TMath::Power(c*d,4.)) + TMath::Power(a*b*c,2.)*(d/TMath::Power(c*d,4.));
81190958 4713
cb690a12 4714 return TMath::Sqrt(error2);
4715}
4716
4717//_____________________________________________________________________________
4718TH1* AliAnalysisMuMu::ComputeEquNofMB(const char* what,const char* quantity,const char* flavour,Bool_t printout)
4719{
81190958 4720
cb690a12 4721 AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
4722 TObjArray* dNchdEtas = binning->CreateBinObjArray();
81190958 4723
cb690a12 4724 Double_t* binArray = binning->CreateBinArray();
81190958 4725
cb690a12 4726 TIter next(dNchdEtas);
4727 AliAnalysisMuMuBinning::Range* r;
81190958 4728
cb690a12 4729 TH1* hFNorm = ComputeDiffFnormFromHistos(what,quantity,flavour,kFALSE);
81190958 4730
cb690a12 4731 TH1* hNMB = new TH1F("hNofEqMB","Equivalent number of MB triggers vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",dNchdEtas->GetEntries(),binArray);
4732
4733 Int_t bin(0);
4734 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
4735 {
4736
4737 TH1* hCMUL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMUL7-B-NOPF-MUON/V0A/%s",
4738 Form("EventsIn%s",r->AsString().Data())));
4739 if ( !hCMUL )
4740 {
4741 AliError(Form("No event histo in bin %s found for CMUL7-B-NOPF-MUON",r->AsString().Data()));
4742 return 0x0;
4743 }
4744
4745 Double_t NMB = hCMUL->GetBinContent(1)*hFNorm->GetBinContent(++bin);
4746 Double_t NMBError = TMath::Sqrt(TMath::Power(hCMUL->GetBinContent(1)*hFNorm->GetBinError(bin),2.) + TMath::Power(TMath::Sqrt(hCMUL->GetBinContent(1))*hFNorm->GetBinContent(bin),2));
4747
4748 if ( printout ) std::cout << r->AsString().Data() << " : " << NMB << " +- " << NMBError << std::endl;
4749
4750 hNMB->SetBinContent(bin,NMB);
4751 hNMB->SetBinError(bin,NMBError);
4752 }
4753
4754 delete dNchdEtas;
4755 delete[] binArray;
4756
4757 return hNMB;
81190958 4758}
4759
cb690a12 4760
a58729a5 4761//_____________________________________________________________________________
1afce1ce 4762AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
2f331ac9 4763{
1afce1ce 4764 /// Correct one spectra
a58729a5 4765
1afce1ce 4766 if (!SIM())
4767 {
4768 AliError("Cannot compute corrected yield without associated MC file !");
4769 return 0x0;
4770 }
4771
cb690a12 4772 const char* accEffSubResultName="PSICOUNT:1";
a58729a5 4773
1afce1ce 4774 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
4775 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
a58729a5 4776
1afce1ce 4777 if ( !realSpectra )
4778 {
4779 AliError("could not get real spectra");
4780 return 0x0;
4781 }
a58729a5 4782
1afce1ce 4783 if ( !simSpectra)
a58729a5 4784 {
1afce1ce 4785 AliError("could not get sim spectra");
a58729a5 4786 return 0x0;
4787 }
4788
1afce1ce 4789 realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
4790
4791 Update();
4792
4793 return realSpectra;
4794}
4795
4796//_____________________________________________________________________________
4797AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
4798{
4799 if (!SIM())
4800 {
4801 AliError("Cannot compute corrected yield without associated MC file !");
4802 return 0x0;
4803 }
4804
cb690a12 4805 const char* accEffSubResultName="PSICOUNT:1";
1afce1ce 4806
4807 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
4808
4809 if ( !realSpectra )
4810 {
4811 AliError("could not get real spectra");
4812 return 0x0;
4813 }
4814
4815 if (!realSpectra->HasValue("CoffNofJpsi"))
4816 {
4817 if (!CorrectSpectra(type,flavour))
4818 {
4819 AliError("Could not get corrected spectra");
4820 return 0x0;
4821 }
4822 }
4823
4824 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
4825
4826 if ( !simSpectra)
2f331ac9 4827 {
a58729a5 4828 AliErrorClass("could not get sim spectra");
4829 return 0x0;
2f331ac9 4830 }
a58729a5 4831
1afce1ce 4832 Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
4833 Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
4834 Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
a58729a5 4835
1afce1ce 4836 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
4837 TObjArray* bins = binning->CreateBinObjArray();
4838 TIter nextBin(bins);
4839 AliAnalysisMuMuBinning::Range* bin;
4840 Int_t i(0);
81190958 4841 AliAnalysisMuMuJpsiResult* r;
a58729a5 4842
1afce1ce 4843 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
4844 {
81190958 4845 r = static_cast<AliAnalysisMuMuJpsiResult*>(realSpectra->BinContentArray()->At(i));
1afce1ce 4846
4847 StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
4848
81190958 4849 AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
1afce1ce 4850
4851 Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
4852 Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4853 nofCINT7,TMath::Sqrt(nofCINT7),
4854 nofCMUL7,TMath::Sqrt(nofCMUL7));
4855
81190958 4856 r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
1afce1ce 4857 nofCINT7,TMath::Sqrt(nofCINT7)));
4858
4859 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
4860
4861 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
4862 mbeq,mbeqError);
4863
4864 r->Set("YJpsi",yield,yieldError);
4865
4866 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
4867 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
4868
4869 ++i;
4870 }
a58729a5 4871
1afce1ce 4872 delete bins;
4873
4874 Update();
4875
4876 return realSpectra;
2f331ac9 4877}
4878
4879//_____________________________________________________________________________
cb690a12 4880AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* type, const char* direction)
2f331ac9 4881{
a58729a5 4882 /// Compute the RAB...
cb690a12 4883
4884 if (!SIM()) return 0x0;
4885
a58729a5 4886 Double_t rapidityShift = 0.465;// 0.5*TMath::Log(208.0/82.0);
4887 const Double_t sqrts=5.023;
4888 const Double_t ymax=TMath::Log(sqrts*1000.0/3.096916);
4889 const Double_t tab = 0.093e-6; // nb^-1
4890 const Double_t tabError = 0.0035E-6; // nb^-1
cb690a12 4891 const char* accEffSubResultName="PSICOUNT:1";
2f331ac9 4892
a58729a5 4893 TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
4894 ydist.SetParameter(0,1.);
4895
4896 //Normalization to the values presented by Zaida and Rosana on January 11th 2013 https://indico.cern.ch/conferenceDisplay.py?confId=224985 slide 22
4897 // Normalization is done in the rapidity range 2.75<y<3.25 where Rosanas values is 230.8+212.1
4898 Double_t y1_norma= 2.75/ymax;
4899 Double_t y2_norma= 3.25/ymax;
4900 Double_t normalization = 0.25*(230.8+212.1)/ydist.Integral(y1_norma, y2_norma);
4901 ydist.SetParameter(0,normalization);
4902// AliInfoClass(Form("ymax=%e normalization=%f",ymax,ydist.Integral(y1_norma, y2_norma)));
2f331ac9 4903
cb690a12 4904 AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
4905 AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(SIM()->OC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
2f331ac9 4906
a58729a5 4907 if ( !realSpectra )
4908 {
4909 AliErrorClass("could not get real spectra");
4910 return 0x0;
4911 }
2f331ac9 4912
a58729a5 4913 if ( !simSpectra)
4914 {
4915 AliErrorClass("could not get sim spectra");
4916 return 0x0;
4917 }
4918
4919 AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
4920 corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
4921
cb690a12 4922 Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
4923 Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
4924 Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
a58729a5 4925
4926 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
4927 TObjArray* bins = binning->CreateBinObjArray();
4928 TIter nextBin(bins);
4929 AliAnalysisMuMuBinning::Range* bin;
4930 Int_t i(0);
81190958 4931 AliAnalysisMuMuJpsiResult* r;
a58729a5 4932
4933 Int_t n = bins->GetLast();
4934
4935 TObjArray finalBins(n+1);
4936 finalBins.SetOwner(kTRUE);
2f331ac9 4937
a58729a5 4938 TObjArray finalResults(n+1);
4939 finalResults.SetOwner(kFALSE);
4940
4941 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
2f331ac9 4942 {
a58729a5 4943 Double_t ylowlab = bin->Xmin();
4944 Double_t yhighlab = bin->Xmax();
4945
4946 Double_t ylowcms, yhighcms;
4947 Double_t ylownorm, yhighnorm;
2f331ac9 4948
5376e016 4949 if ( bin->IsIntegrated() )
a58729a5 4950 {
4951 ylowlab = -4;
4952 yhighlab = -2.5;
4953 }
2f331ac9 4954
a58729a5 4955 if ( strcmp(direction,"pPb")==0 )
2f331ac9 4956 {
a58729a5 4957 ylowcms = TMath::Abs(yhighlab) - rapidityShift;
4958 yhighcms = TMath::Abs(ylowlab) - rapidityShift;
4959 ylownorm = ylowcms/ymax;
4960 yhighnorm = yhighcms/ymax;
2f331ac9 4961 }
4962 else
4963 {
a58729a5 4964 ylowcms = ylowlab - rapidityShift;
4965 yhighcms = yhighlab - rapidityShift;
4966 ylownorm = -yhighcms/ymax;
4967 yhighnorm = -ylowcms/ymax;
2f331ac9 4968 }
4969
2f331ac9 4970
a58729a5 4971 Double_t brsigmapp = ydist.Integral(ylownorm,yhighnorm);
4972 Double_t brsigmappError = 0.0; // FIXME
4973
4974 AliInfoClass(Form("y range : LAB %f ; %f CMS %f ; %f -> ynorm : %f ; %f -> BR x sigmapp = %f",
4975 ylowlab,yhighlab,ylowcms,yhighcms,ylownorm,yhighnorm,brsigmapp));
4976
81190958 4977 r = static_cast<AliAnalysisMuMuJpsiResult*>(corrSpectra->BinContentArray()->At(i)->Clone());
a58729a5 4978
81190958 4979 AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
a58729a5 4980
4981 Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
4982 Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4983 nofCINT7,TMath::Sqrt(nofCINT7),
4984 nofCMUL7,TMath::Sqrt(nofCMUL7));
4985
81190958 4986 r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
a58729a5 4987 nofCINT7,TMath::Sqrt(nofCINT7)));
4988
4989 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
4990
4991 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
4992 mbeq,mbeqError);
4993
4994 r->Set(Form("Y%sJpsi",direction),yield,yieldError);
4995
4996 Double_t raa = yield/(tab*brsigmapp);
4997 Double_t raaError = AliAnalysisMuMuResult::ErrorABC(yield,yieldError,
4998 tab,tabError,
4999 brsigmapp,brsigmappError);
5000 r->Set(Form("R%sJpsi",direction),raa,raaError);
5001
5002 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
5003 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
5004
5376e016 5005 AliAnalysisMuMuBinning::Range* bincm = new AliAnalysisMuMuBinning::Range(bin->What(),bin->Quantity(),ylowcms,yhighcms);
a58729a5 5006
5007 r->SetBin(*bincm);
5008
5009 finalBins.Add(bincm);
5010 finalResults.Add(r);
5011
5012 ++i;
2f331ac9 5013 }
5014
a58729a5 5015 delete bins;
2f331ac9 5016
a58729a5 5017 AliAnalysisMuMuSpectra* spectra = new AliAnalysisMuMuSpectra(type,direction);
5018
5019 for ( i = 0; i <= n; ++i )
5020 {
5021 Int_t j(i);
5022 if ( strcmp(direction,"pPb")==0 )
5023 {
5024 j = n-i;
5025 }
5026
81190958 5027 r = static_cast<AliAnalysisMuMuJpsiResult*>(finalResults.At(j));
a58729a5 5028
5029 bin = static_cast<AliAnalysisMuMuBinning::Range*>(finalBins.At(j));
5030
5031 spectra->AdoptResult(*bin,r);
5032 }
5033
5034
5035 delete corrSpectra;
2f331ac9 5036
a58729a5 5037 return spectra;
2f331ac9 5038}
5039
5040//_____________________________________________________________________________
cb690a12 5041void AliAnalysisMuMu::SetConfig(const AliAnalysisMuMuConfig& config)
2f331ac9 5042{
cb690a12 5043 /// (re)set the config
5044 delete fConfig;
5045 fConfig = new AliAnalysisMuMuConfig(config);
2f331ac9 5046}
a58729a5 5047