]> 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
Fix compilation warnings (Laurent)
[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"
20
a58729a5 21#include "AliAnalysisMuMuBinning.h"
22#include "AliAnalysisMuMuResult.h"
23#include "AliAnalysisMuMuSpectra.h"
2f331ac9 24#include "AliAnalysisTriggerScalers.h"
25#include "AliCounterCollection.h"
26#include "AliHistogramCollection.h"
27#include "AliLog.h"
a58729a5 28#include "AliMergeableCollection.h"
2f331ac9 29#include "Riostream.h"
30#include "TArrayL64.h"
a58729a5 31#include "TASImage.h"
32#include "TAxis.h"
33#include "TCanvas.h"
34#include "TColor.h"
35#include "TF1.h"
36#include "TFile.h"
37#include "TGraphErrors.h"
38#include "TGrid.h"
39#include "TH1.h"
40#include "TH2.h"
41#include "TKey.h"
42#include "TLegend.h"
43#include "TLegendEntry.h"
44#include "TLine.h"
45#include "TList.h"
46#include "TMap.h"
47#include "TMath.h"
48#include "TObjArray.h"
49#include "TObjString.h"
50#include "TParameter.h"
51#include "TPaveText.h"
52#include "TROOT.h"
1afce1ce 53#include "TStopwatch.h"
a58729a5 54#include "TStyle.h"
55#include "TSystem.h"
56#include <cassert>
57#include <map>
58#include <set>
59#include <string>
2f331ac9 60
a58729a5 61ClassImp(AliAnalysisMuMu)
2f331ac9 62
a58729a5 63TString AliAnalysisMuMu::fgOCDBPath("raw://");
2f331ac9 64
1afce1ce 65TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-B-NOPF-MUON");
66
67//,CMUL7-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-MUON,CMUL8-S-NOPF-MUON,CMUL7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-MUON,CPBI1MUL-B-NOPF-MUON,CMULLO-B-NOPF-MUON");
68
69TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON");
2f331ac9 70
1afce1ce 71//,CMSL7-S-NOPF-ALLNOTRD,CMSL8-S-NOPF-MUON,CMSL8-S-NOPF-ALLNOTRD,CMSL7-B-NOPF-MUON,CMUS1-B-NOPF-MUON,CMUS7-B-NOPF-MUON,CMSNGL-B-NOPF-MUON");
2f331ac9 72
1afce1ce 73TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD");
2f331ac9 74
1afce1ce 75//,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
76
77TString AliAnalysisMuMu::fgDefaultEventSelectionList("PSALL"); // for real data, for simulation see AliAnalysisMuMu ctor
2f331ac9 78
a58729a5 79TString AliAnalysisMuMu::fgDefaultPairSelectionList("pMATCHLOWRABSBOTH");
2f331ac9 80
a58729a5 81TString AliAnalysisMuMu::fgDefaultCentralitySelectionList("PP");
2f331ac9 82
1afce1ce 83//TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOW:2,PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
84TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
85
86TString AliAnalysisMuMu::fgDefaultEventSelectionForSimulations("ALL");
87TString AliAnalysisMuMu::fgDefaultDimuonTriggerForSimulations("CMULLO-B-NOPF-MUON");
88
a58729a5 89Bool_t AliAnalysisMuMu::fgIsCompactGraphs(kFALSE);
2f331ac9 90
1afce1ce 91//_____________________________________________________________________________
92TString First(const TString s)
93{
94 TString rv;
95
96 TObjArray* tokens = s.Tokenize(",");
97
98 if (!tokens) return rv;
99
100 rv = static_cast<TObjString*>(tokens->First())->String();
101
102 delete tokens;
103
104 return rv;
105}
106
2f331ac9 107//_____________________________________________________________________________
a58729a5 108TString FindTrigger(const AliMergeableCollection& mc,
2f331ac9 109 const char* base,
110 const char* selection,
111 const char* paircut,
112 const char* centrality)
113{
114 /// find the trigger containing the MinvPt histograms
115
116 std::vector<std::string> trigger2test;
117
118 // trigger2test.push_back(Form("%s5-B-NOPF-ALLNOTRD",base));
119 // trigger2test.push_back(Form("%s1-B-NOPF-ALLNOTRD",base));
120 // trigger2test.push_back(Form("%s1B-ABCE-NOPF-MUON",base));
121 if ( TString(base).Contains("||") || TString(base).Contains("-") )
122 {
123 trigger2test.push_back(base);
124 }
125 else
126 {
127 trigger2test.push_back(Form("%s-B-NOPF-ALLNOTRD",base));
128 trigger2test.push_back(Form("%s-B-NOPF-MUON",base));
129 trigger2test.push_back(Form("%s-S-NOPF-ALLNOTRD",base));
130 trigger2test.push_back(Form("%s-S-NOPF-MUON",base));
131 }
132 trigger2test.push_back("ANY");
133
a58729a5 134 for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
2f331ac9 135 {
136 std::string trigger = trigger2test[i];
137
a58729a5 138 if ( mc.GetObject(Form("/%s/%s/%s/%s",selection,trigger.c_str(),centrality,paircut),"MinvUS") ||
139 mc.GetObject(Form("/%s/%s/%s/%s",selection,trigger.c_str(),centrality,paircut),"MinvUSPt")
140 )
2f331ac9 141 {
142 return trigger.c_str();
143 }
144 }
145
a58729a5 146 // AliWarningGeneral("FindTrigger",Form("DID NOT FIND TRIGGER base=%s selection=%s paircut=%s centrality=%s",
147 // base,selection,paircut,centrality));
148 // for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
149 // {
150 // AliWarningGeneral("FindTrigger",Form("tested trigger = %s",trigger2test[i].c_str()));
151 // }
2f331ac9 152 return "";
153}
154
2f331ac9 155//_____________________________________________________________________________
1afce1ce 156AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName) : TObject(),
2f331ac9 157fFilename(filename),
2f331ac9 158fCounterCollection(0x0),
2754fd07 159fDimuonTriggers(fgDefaultDimuonTriggers),
160fMuonTriggers(fgDefaultMuonTriggers),
161fMinbiasTriggers(fgDefaultMinbiasTriggers),
162fEventSelectionList(fgDefaultEventSelectionList),
a58729a5 163fPairSelectionList(fgDefaultPairSelectionList),
164fCentralitySelectionList(fgDefaultCentralitySelectionList),
1afce1ce 165fFitTypeList(fgDefaultFitTypeList),
a58729a5 166fBinning(0x0),
167fMergeableCollection(0x0),
168fRunNumbers(),
1afce1ce 169fCorrectionPerRun(0x0),
170fAssociatedSimulation(0x0)
2f331ac9 171{
172 // ctor
173
a58729a5 174 GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
1afce1ce 175
176 if (IsSimulation())
177 {
178 SetEventSelectionList("ALL");
179 SetDimuonTriggerList("CMULLO-B-NOPF-MUON");
180 SetFitTypeList("PSI1:1,COUNTJPSI:1");
181 }
182
183 if ( strlen(associatedSimFileName) )
184 {
185 fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
186 }
2f331ac9 187}
188
189//_____________________________________________________________________________
190AliAnalysisMuMu::~AliAnalysisMuMu()
191{
192 // dtor
2f331ac9 193 delete fCounterCollection;
a58729a5 194 delete fBinning;
195 delete fMergeableCollection;
196 delete fCorrectionPerRun;
1afce1ce 197 delete fAssociatedSimulation;
2f331ac9 198}
199
200//_____________________________________________________________________________
2754fd07 201void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
202 ULong64_t* totalNmb,
203 ULong64_t* totalNmsl,
204 ULong64_t* totalNmul)
2f331ac9 205{
206 // Report of some basic numbers, like number of MB and MUON triggers,
207 // both before and after physics selection, and comparison with
208 // the total number of such triggers (taken from the OCDB scalers)
209 // if requested.
210 //
211 // filename is assumed to be a root filecontaining a list containing
212 // an AliCounterCollection (or directly an AliCounterCollection)
213 //
214 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
215 //
216
a58729a5 217 if (!fMergeableCollection || !fCounterCollection) return;
2f331ac9 218
219 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
220 TIter nextRun(runs);
221
222 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
223 TIter nextTrigger(triggers);
224
225 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
226
2754fd07 227 Bool_t doPS = (events->FindObject("PSALL") != 0x0);
228
2f331ac9 229 TObjString* srun;
230 TObjString* strigger;
231
2754fd07 232 ULong64_t localNmb(0);
233 ULong64_t localNmsl(0);
234 ULong64_t localNmul(0);
235
236 if ( totalNmb) *totalNmb = 0;
237 if ( totalNmsl) *totalNmsl = 0;
238 if ( totalNmul ) *totalNmul = 0;
2f331ac9 239
240 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
241 {
242 std::cout << Form("RUN %09d ",srun->String().Atoi());
243
244 TString details;
245 ULong64_t nmb(0);
246 ULong64_t nmsl(0);
247 ULong64_t nmul(0);
248
249 nextTrigger.Reset();
250
2754fd07 251 Int_t nofPS(0);
252
2f331ac9 253 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
254 {
2754fd07 255
256 if ( !fgDefaultMinbiasTriggers.Contains(strigger->String().Data()) &&
257 !fgDefaultMuonTriggers.Contains(strigger->String().Data()) &&
258 !fgDefaultDimuonTriggers.Contains(strigger->String().Data()) ) continue;
259
81dfe194 260 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
261 strigger->String().Data(),"ALL",srun->String().Atoi())));
2f331ac9 262
2754fd07 263 details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
2f331ac9 264
265
81dfe194 266 ULong64_t nps = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
267 strigger->String().Data(),"PSALL",srun->String().Atoi())));
2754fd07 268
269 if ( doPS )
2f331ac9 270 {
2754fd07 271 details += TString::Format(" PS %5.1f %%",nps*100.0/n);
2f331ac9 272 }
273
2754fd07 274 if (nps)
275 {
276 ++nofPS;
277 }
278
2f331ac9 279 if ( fMinbiasTriggers.Contains(strigger->String()) )
280 {
281 nmb += n;
2754fd07 282 if ( totalNmb) (*totalNmb) += n;
283 localNmb += n;
2f331ac9 284 }
285 else if ( fMuonTriggers.Contains(strigger->String()) )
286 {
287 nmsl += n;
2754fd07 288 if ( totalNmsl) (*totalNmsl) += n;
289 localNmsl += n;
2f331ac9 290 }
291 else if ( fDimuonTriggers.Contains(strigger->String()) )
292 {
293 nmul += n;
2754fd07 294 if ( totalNmul ) (*totalNmul) += n;
295 localNmul += n;
2f331ac9 296 }
297 }
298
2754fd07 299 std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
300 nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
2f331ac9 301
302 if ( detailTriggers )
303 {
304 std::cout << details.Data();
305 }
306 std::cout << std::endl;
307 }
308
2754fd07 309 if ( !totalNmul && !totalNmsl && !totalNmb )
310 {
311 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
312 localNmb,localNmsl,localNmul) << std::endl;
313 }
2f331ac9 314
2754fd07 315 delete runs;
2f331ac9 316 delete triggers;
317 delete events;
318}
319
2754fd07 320//_____________________________________________________________________________
321void AliAnalysisMuMu::BasicCountsEvolution(const char* filelist, Bool_t detailTriggers)
322{
323 // Report of some basic numbers, like number of MB and MUON triggers,
324 // both before and after physics selection, and comparison with
325 // the total number of such triggers (taken from the OCDB scalers)
326 // if requested.
327 //
328 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
329 //
330 // To change the list of (single muon, dimuon, MB) triggers, use
331 // the SetDefault*TriggerList methods prior to call this one
332 //
333
334 TObjArray* files = ReadFileList(filelist);
335
336 if (!files || files->IsEmpty() ) return;
337
338 TIter next(files);
339 TObjString* str;
340
341 ULong64_t totalNmb(0);
342 ULong64_t totalNmsl(0);
343 ULong64_t totalNmul(0);
344
345 while ( ( str = static_cast<TObjString*>(next()) ) )
346 {
347 AliAnalysisMuMu m(str->String().Data());
348
349 ULong64_t nmb(0);
350 ULong64_t nmsl(0);
351 ULong64_t nmul(0);
352
353 m.BasicCounts(detailTriggers,&nmb,&nmsl,&nmul);
354
355 totalNmb += nmb;
356 totalNmsl += nmsl;
357 totalNmul += nmul;
358 }
359
360 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
361 totalNmb,totalNmsl,totalNmul) << std::endl;
362
363}
364
a58729a5 365//_____________________________________________________________________________
366void AliAnalysisMuMu::CentralityCheck(const char* filelist)
367{
368 // Check if we get correctly filled centrality
369
370 TObjArray* files = ReadFileList(filelist);
371
372 if (!files || files->IsEmpty() ) return;
373
374 TIter next(files);
375 TObjString* str;
376
377 while ( ( str = static_cast<TObjString*>(next()) ) )
378 {
379 AliMergeableCollection* mc(0x0);
380 AliCounterCollection* cc(0x0);
381 AliAnalysisMuMuBinning* bin(0x0);
382 std::set<int> runnumbers;
383
384 if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
385
386 int run = RunNumberFromFileName(str->String().Data());
387
388 TH1* h = mc->Histo("/ALL/CPBI1MUL-B-NOPF-MUON/Centrality");
389
390 float percent(0);
391
392 if (h)
393 {
394 percent = 100*h->Integral(1,1)/h->Integral();
395 }
396
397 std::cout << Form("RUN %09d PERCENT %7.2f",run,percent) << std::endl;
398
399 delete mc;
400 }
401
402 gROOT->CloseFiles();
403
404}
405
1afce1ce 406//_____________________________________________________________________________
407void AliAnalysisMuMu::CleanAllSpectra()
408{
409 /// Delete all the spectra we may have
410
411 MC()->RemoveByType("AliAnalysisMuMuSpectra");
412 Update();
413}
a58729a5 414
415//_____________________________________________________________________________
416void AliAnalysisMuMu::Compact(TGraph& g)
417{
418 /// Compact (i.e. get the equivalent of 1 bin = 1 run number for an histogram)
419 /// the graph.
420 /// Only works if the x content of this graph represents run numbers. Otherwise
421 /// result is unpredictable ;-)
422
423 if ( !g.GetN() ) return;
424
425 TGraph* newgraph = static_cast<TGraph*>(g.Clone());
426 Double_t x,xerr,y,yerr;
427 TGraphErrors* ge = dynamic_cast<TGraphErrors*>(newgraph);
428
429 TAxis* axis = g.GetXaxis();
430
431 for ( Int_t i = 0; i < newgraph->GetN(); ++i )
432 {
433 g.GetPoint(i,x,y);
434 if (ge)
435 {
436 xerr = ge->GetErrorX(i);
437 yerr = ge->GetErrorY(i);
438 }
439
440 g.SetPoint(i,i+0.5,y);
441 if (ge)
442 {
443 static_cast<TGraphErrors&>(g).SetPointError(i,0.5,yerr);
444 }
445
446 axis->SetBinLabel(i,Form("%d",TMath::Nint(x)));
447 }
448
449
450}
451
452
2f331ac9 453//_____________________________________________________________________________
454TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
455 const char* backgroundresults)
456{
457 TFile* fjpsi = FileOpen(jpsiresults);
458 TFile* fbck = FileOpen(backgroundresults);
459
460 if (!fjpsi || !fbck) return 0x0;
461
462 TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
463
464 std::vector<std::string> checks;
465
466 checks.push_back("muminus-CMUU7-B-NOPF-ALLNOTRD");
467 checks.push_back("muplus-CMUU7-B-NOPF-ALLNOTRD");
468 checks.push_back("muminus-CMUSH7-B-NOPF-MUON");
469 checks.push_back("muplus-CMUSH7-B-NOPF-MUON");
470
471 if (!gjpsi) return 0x0;
472
473 TObjArray* a = new TObjArray;
474 a->SetOwner(kTRUE);
475
476 for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
477 {
478
479 TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
480
481 if (!gback) continue;
482
483 if ( gjpsi->GetN() != gback->GetN() )
484 {
485 AliErrorClass("graphs have different number of points !");
486 continue;
487 }
488
489 TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
490
491 for ( int i = 0; i < gjpsi->GetN(); ++i )
492 {
493 double r1,r2,y1,y2;
494
495 gjpsi->GetPoint(i,r1,y1);
496 gback->GetPoint(i,r2,y2);
497
498 if ( r1 != r2 )
499 {
500 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
501 continue;
502 }
503
504 g->SetPoint(i,y2,y1);
505 // g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
506 }
507
508 g->SetMarkerStyle(25+j);
509 g->SetMarkerSize(1.2);
510 if (j==0)
511 {
512 g->Draw("ap");
513 }
514 else
515 {
516 g->Draw("p");
517 }
518 g->SetLineColor(j+1);
519 g->SetMarkerColor(j+1);
520 g->SetName(checks[j].c_str());
521 a->AddLast(g);
522 }
523
524 return a;
525}
526
527//_____________________________________________________________________________
528TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
529 const char* simjpsiresults)
530{
531 TFile* freal = FileOpen(realjpsiresults);
532 TFile* fsim = FileOpen(simjpsiresults);
533
534 if (!freal || !fsim) return 0x0;
535
536 TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
537 TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
538
539 TObjArray* a = new TObjArray;
540 a->SetOwner(kTRUE);
541
542 if ( greal->GetN() != gsim->GetN() )
543 {
544 AliErrorClass("graphs have different number of points !");
545 return 0x0;
546 }
547
548 TGraphErrors* g = new TGraphErrors(greal->GetN());
549 TGraphErrors* gratio = new TGraphErrors(greal->GetN());
550
551 for ( int i = 0; i < greal->GetN(); ++i )
552 {
553 double r1,r2,y1,y2;
554
555 greal->GetPoint(i,r1,y1);
556 gsim->GetPoint(i,r2,y2);
557
558 if ( r1 != r2 )
559 {
560 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
561 continue;
562 }
563
564 double ratio(0.0);
565
566 if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
567 {
568 g->SetPoint(i,0,0);
569 g->SetPointError(i,0,0);
570 }
571 else
572 {
573 g->SetPoint(i,y2,y1);
574 g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
575 ratio = y2/y1;
576 }
577 gratio->SetPoint(i,r1,ratio);
578 }
579
580 g->SetMarkerStyle(25);
581 g->SetMarkerSize(1.2);
582
583 new TCanvas;
584
585 g->Draw("ap");
586
587 g->SetLineColor(1);
588 g->SetMarkerColor(1);
589 g->SetName("jpsipercmuurealvssim");
590
591 new TCanvas;
592
593 greal->Draw("alp");
594 gsim->SetLineColor(4);
595
596 gsim->Draw("lp");
597
598 new TCanvas;
599 gratio->Draw("alp");
600
601 return g;
602}
603
604//_____________________________________________________________________________
605TObjArray* AliAnalysisMuMu::ComputeBackgroundEvolution(const char* filelist,
a58729a5 606 const char* triggerList,
607 Double_t ptmin,
2f331ac9 608 const char* outputFile,
609 const char* outputMode)
610{
611 // triggerList is a list of complete trigger names, separated by space
612 // of the triggers to consider : only the first one found in the list
613 // is used for each run.
614
615 TObjArray* files = ReadFileList(filelist);
616
617 if (!files || files->IsEmpty() ) return 0x0;
618
619 TIter next(files);
620 TObjString* str;
621
a58729a5 622 const char* ps = "PSALL";
2f331ac9 623 const char* centrality = "PP";
624 const char* ts1 = "sMATCHLOWRABS";
625 const char* ts2 = "sMATCHLOWRABSDCA";
626
627 std::map<std::string, std::vector<float> > runs;
628 std::map<std::string, std::vector<float> > errruns;
629 std::map<std::string, std::vector<float> > yplus,erryplus;
630 std::map<std::string, std::vector<float> > yminus,erryminus;
631
2754fd07 632 TObjArray* triggers = TString(triggerList).Tokenize(",");
2f331ac9 633
634 TObjArray* a = new TObjArray;
635 a->SetOwner(kTRUE);
636
637 Bool_t bothSigns(kFALSE);
638
639 while ( ( str = static_cast<TObjString*>(next()) ) )
640 {
641 AliInfoClass(str->String().Data());
642
a58729a5 643 AliMergeableCollection* mc(0x0);
2f331ac9 644 AliCounterCollection* cc(0x0);
a58729a5 645 AliAnalysisMuMuBinning* bin(0x0);
646 std::set<int> runnumbers;
2f331ac9 647
a58729a5 648 if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
2f331ac9 649
a58729a5 650 TIter nextObject(mc->CreateIterator());
651 TObject* o;
2f331ac9 652 int nplus(0), nminus(0);
653
a58729a5 654 while ( ( o = nextObject() ) )
2f331ac9 655 {
a58729a5 656 if ( o->InheritsFrom("TH1") )
657 {
658 continue;
659 }
660
661 TH1* h = static_cast<TH1*>(o);
662
2f331ac9 663 if ( TString(h->GetName()).EndsWith("Plus") )
664 {
665 nplus++;
666 }
667 if ( TString(h->GetName()).EndsWith("Minus") )
668 {
669 nminus++;
670 }
671 }
672
673 if (nminus==nplus && nplus>0 )
674 {
675 bothSigns = kTRUE;
676 }
677
678 AliInfoClass(Form("Both signs = %d",bothSigns));
679
680 TIter nextTrigger(triggers);
681 TObjString* trigger;
682 TH1* h1p(0x0);
683 TH1* h1m(0x0);
684 TH1* h2p(0x0);
685 TH1* h2m(0x0);
686 TString format;
687
688 while ( ( trigger = static_cast<TObjString*>(nextTrigger()) ) )
689 {
690 if (bothSigns)
691 {
692 format = "/%s/%s/%s/%s/PtEtaMuPlus:py";
693 }
694 else
695 {
696 format = "/%s/%s/%s/%s/PtEtaMu:py";
697 }
698
a58729a5 699 TString hname(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
700
701 h1p = mc->Histo(hname.Data());
2f331ac9 702
703 if (!h1p)
704 {
a58729a5 705 AliInfoClass(Form("Could not get %s",hname.Data()));
2f331ac9 706 continue;
707 }
708
709 AliInfoClass(Form("Will use trigger %s",trigger->String().Data()));
710
a58729a5 711 h2p = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
2f331ac9 712
713 if ( bothSigns )
714 {
715 format.ReplaceAll("Plus","Minus");
a58729a5 716 h1m = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
717 h2m = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
2f331ac9 718 }
719 else
720 {
721 h2m=h2p;
722 h1m=h1p;
723 }
724
725 if (h1m && h2m && h1p && h2p)
726 {
a58729a5 727 Int_t bin1 = h1m->GetXaxis()->FindBin(ptmin);
728 Int_t bin2 = h1m->GetXaxis()->GetNbins();
729
2f331ac9 730 runs[trigger->String().Data()].push_back(RunNumberFromFileName(str->String().Data()));
731 errruns[trigger->String().Data()].push_back(0.5);
a58729a5 732
733 double e1,e2;
734 double v1 = h2m->IntegralAndError(bin1,bin2,e1);
735 double v2 = h1m->IntegralAndError(bin1,bin2,e2);
736 double value = 100*(1.0-v1/v2);
737 e1/=v1;
738 e2/=v2;
2f331ac9 739 yminus[trigger->String().Data()].push_back(value);
a58729a5 740// double e1 = 1.0/TMath::Sqrt(h1m->GetEntries());
741// double e2 = 1.0/TMath::Sqrt(h2m->GetEntries());
2f331ac9 742 erryminus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
743
a58729a5 744 v1=h2p->IntegralAndError(bin1,bin2,e1);
745 v2=h1p->IntegralAndError(bin1,bin2,e1);
746 value = 100*(1.0-v1/v2);
747 e1/=v1;
748 e2/=v2;
2f331ac9 749 yplus[trigger->String().Data()].push_back(value);
a58729a5 750// e1 = 1.0/TMath::Sqrt(h1p->GetEntries());
751// e2 = 1.0/TMath::Sqrt(h2p->GetEntries());
2f331ac9 752 erryplus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
753 }
754 else
755 {
756 std::cout << Form("Error : h1m %p h2m %p h1p %p h2p %p",h1m,h2m,h1p,h2p) << std::endl;
757 }
758 }
759
a58729a5 760 delete mc;
2f331ac9 761 delete cc;
762 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(str->String().Data()));
763 delete f;
764 }
765
766 delete triggers;
767
768 TFile* f = new TFile(outputFile,outputMode);
769
770 std::map<std::string, std::vector<float> >::const_iterator it;
771
772 for ( it = runs.begin(); it != runs.end(); ++it )
773 {
774 std::string triggerName = it->first;
775
776 TGraphErrors* gp = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yplus[triggerName][0],&errruns[triggerName][0],&erryplus[triggerName][0]);
777 TGraphErrors* gm(0x0);
778
779 if ( bothSigns )
780 {
781 gm = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yminus[triggerName][0],&errruns[triggerName][0],&erryminus[triggerName][0]);
782 }
783
784 if ( bothSigns )
785 {
786 gp->Write(Form("muplus_%s",triggerName.c_str()),TObject::kOverwrite);
787 gm->Write(Form("muminus_%s",triggerName.c_str()),TObject::kOverwrite);
788 }
789 else
790 {
791 gp->Write(Form("mu_%s",triggerName.c_str()),TObject::kOverwrite);
792 }
793
794 }
795
796 delete f;
797
798 return a;
799}
800
2754fd07 801//_____________________________________________________________________________
802TMap*
803AliAnalysisMuMu::ComputeJpsiEvolution(const char* filelist, const char* triggerList,
a58729a5 804 const char* outputFile)
2754fd07 805{
806 /// Compute some jpsi information for a list of files / trigger combinations
807
808 TObjArray* files = ReadFileList(filelist);
809
810 if (!files || files->IsEmpty() ) return 0x0;
811
812 TMap results; // one TObjString->TObjArray per file
813 results.SetOwnerKeyValue(kTRUE,kTRUE);
814
815 TIter nextFile(files);
816 TObjString* str;
a58729a5 817 TString fitType;
2754fd07 818
a58729a5 819// while ( ( str = static_cast<TObjString*>(nextFile()) ) )
820// {
821// std::cout << str->String().Data() << std::endl;
822//
823// AliAnalysisMuMu m(str->String().Data());
824//
825// m.SetDimuonTriggerList(triggerList);
826//
827// TMap* map = m.Jpsi();
828//
829// if (!map)
830// {
831// AliWarningClass(Form("Got no jpsi for %s",str->String().Data()));
832// }
833//
834// results.Add(new TObjString(str->String()), map);
835// }
836//
837// if (!results.GetSize()) return 0x0;
2754fd07 838
839 // compute the total over all files
840
841 TMap* total = new TMap;
842 total->SetOwnerKeyValue(kTRUE,kTRUE);
843
844 nextFile.Reset();
845 TObjArray* triggers = TString(triggerList).Tokenize(",");
846
847 TIter nextTrigger(triggers);
848 TObjString* trigger(0x0);
849
850 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
851 {
852 nextFile.Reset();
853
a58729a5 854 TList l;
855 AliAnalysisMuMuResult* ref(0x0);
2754fd07 856
857 while ( ( str = static_cast<TObjString*>(nextFile()) ) )
858 {
a58729a5 859// TObjArray* a = static_cast<TObjArray*>(results.GetValue(str->String().Data()));
2754fd07 860
a58729a5 861 AliInfoClass("FIXME: write the merging of AliAnalysisMuMuResult objects !");
2754fd07 862
a58729a5 863// AliAnalysisMuMuResult* r(0x0);
864//
865// if (a)
866// {
867// r = static_cast<AliAnalysisMuMuResult*>(a->FindObject(trigger->String().Data()));
868//
869// if (r)
870// {
871// if (!ref) ref = r;
872//
873// if ( !hminv )
874// {
875// TH1* htmp = static_cast<TH1*>(r->Minv()->Clone());
876// if (!htmp)
877// {
878// continue;
879// }
880// hminv = htmp;
881// }
882// else
883// {
884// l.Add(r->Minv());
885// }
886//
887// n += r->NofTriggers();
888// ++nruns;
889// }
890// }
2754fd07 891 }
892
a58729a5 893// sum->Merge(&l);
2754fd07 894
895 if (!ref) continue;
896
2754fd07 897
a58729a5 898// AliAnalysisMuMuResult* sum = new AliAnalysisMuMuResult(*hminv,
899// ref->TriggerClass(),
900// ref->EventSelection(),
901// ref->PairSelection(),
902// ref->CentralitySelection(),
903// AliAnalysisMuMuBinning::Range());
904//
905// sum->SetNofTriggers(n);
906//
907// sum->SetNofRuns(nruns);
908//
909// sum->Fit(1);
910//
911// total->Add(new TObjString(trigger->String().Data()),sum);
2754fd07 912
913 }
914
a58729a5 915 AliInfoClass("--------------------------------------");
2754fd07 916 StdoutToAliInfoClass(total->Print(););
a58729a5 917
918 AliInfoClass("---------------------------Going to write file");
919
920 TFile* fout = new TFile(outputFile,"RECREATE");
2754fd07 921
a58729a5 922 results.Write("rbr",TObject::kSingleKey);
2754fd07 923
a58729a5 924 total->Write("total",TObject::kSingleKey);
925
926 fout->Close();
2754fd07 927
928 delete fout;
929
930 AliInfoClass(Form("%d files analyzed",files->GetEntries()));
931
932 return total;
933}
934
2f331ac9 935//_____________________________________________________________________________
936Bool_t AliAnalysisMuMu::DecodeFileName(const char* filename,
937 TString& period,
938 int& esdpass,
939 int& aodtrain,
940 int& runnumber)
941{
942 esdpass=aodtrain=runnumber=-1;
943 period="";
944
945 TString sfile(gSystem->BaseName(filename));
946
947 if (!sfile.BeginsWith("LHC") && !sfile.BeginsWith("SIM") )
948 {
949 std::cerr << Form("filename %s does not start with LHC or SIM",filename) << std::endl;
950 return kFALSE;
951 }
952
953 int year;
954 char p;
955 Bool_t ok(kFALSE);
956
957 if ( sfile.BeginsWith("LHC") )
958 {
959 if ( sfile.Contains("pass") && sfile.Contains("AODMUON") )
960 {
961 int pass;
962 sscanf(sfile.Data(),"LHC%2d%c_pass%d_AODMUON%03d_%09d",&year,&p,&pass,&aodtrain,&runnumber);
963 period = TString::Format("LHC%2d%c",year,p);
964 ok=kTRUE;
965 }
966 else if ( sfile.Contains("pass") && sfile.Contains("_muon_") && sfile.Contains("AOD000") )
967 {
968 // LHC11c_pass2_muon_AOD000_000152087.saf.root
969 sscanf(sfile.Data(),"LHC%2d%c_pass%d_muon_AOD000_%09d",&year,&p,&esdpass,&runnumber);
970 period = TString::Format("LHC%2d%c",year,p);
971 ok=kTRUE;
972 }
973 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AODMUON000") )
974 {
975 // LHC12h_muon_calo_AODMUON000_000190112.saf.root
976 sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AODMUON000_%09d",&year,&p,&runnumber);
977 period = TString::Format("LHC%2d%c",year,p);
978 ok=kTRUE;
979 }
980 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AOD000") )
981 {
982 // LHC12h_muon_calo_AOD000_000190112.saf.root
983 sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AOD000_%09d",&year,&p,&runnumber);
984 period = TString::Format("LHC%2d%c",year,p);
985 ok=kTRUE;
986 }
987 else if ( sfile.Contains("AODMUON" ) )
988 {
989 sscanf(sfile.Data(),"LHC%2d%c_AODMUON%03d_%09d",&year,&p,&aodtrain,&runnumber);
990 period = TString::Format("LHC%2d%c",year,p);
991 ok=kTRUE;
992 }
993 else if ( sfile.Contains("AOD000") )
994 {
995 sscanf(sfile.Data(),"LHC%2d%c_muons_AOD000_%09d",&year,&p,&runnumber);
996 period = TString::Format("LHC%2d%c",year,p);
997 ok=kTRUE;
998 }
999 else if ( sfile.Contains("ESD_OUTER000"))
1000 {
1001 sscanf(sfile.Data(),"LHC%2d%c_cpass1_ESD_OUTER000_%09d",&year,&p,&runnumber);
1002 ok=kTRUE;
1003 }
1004 }
1005 else if ( sfile.BeginsWith("SIM_JPSI3" ) )
1006 {
1007 sscanf(sfile.Data(),"SIM_JPSI3_%09d",&runnumber);
1008 ok = kTRUE;
1009 }
1010 else if ( sfile.BeginsWith("SIM_UPSILON" ) )
1011 {
1012 sscanf(sfile.Data(),"SIM_UPSILON_%09d",&runnumber);
1013 ok = kTRUE;
1014 }
1015 else if ( sfile.BeginsWith("SIM_JPSI" ) )
1016 {
1017 sscanf(sfile.Data(),"SIM_JPSI_LHC%2d%c_%09d",&year,&p,&runnumber);
1018 period = TString::Format("LHC%2d%c",year,p);
1019 ok = kTRUE;
1020 }
1021
1022 if (!ok)
1023 {
1024 std::cerr << Form("Can not decode %s",filename) << std::endl;
1025 return kFALSE;
1026 }
1027
1028 return kTRUE;
1029}
1030
a58729a5 1031//_____________________________________________________________________________
1032void AliAnalysisMuMu::DrawMinv(const char* type,
1033 const char* particle,
1034 const char* trigger,
1035 const char* eventType,
1036 const char* pairCut,
1afce1ce 1037 const char* centrality,
1038 const char* subresultname,
1039 const char* flavour) const
a58729a5 1040{
1041 /// Draw minv spectra for binning of given type
1042
1043 if (!MC() || !BIN()) return;
1044
1afce1ce 1045 TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
a58729a5 1046 if (!bins)
1047 {
1048 AliError(Form("Could not get %s bins",type));
1049 return;
1050 }
1051
1afce1ce 1052 Double_t xmin(-1);
1053 Double_t xmax(-1);
1054
1055 TString sparticle(particle);
1056 if ( sparticle=="PSI" )
1057 {
1058 xmin = 2;
1059 xmax = 6;
1060 }
1061
1062 Int_t nx(1);
1063 Int_t ny(1);
1064
1065 Int_t n = bins->GetEntries();
1066
1067 if ( n == 2 )
1068 {
1069 nx = 2;
1070 }
1071 else if ( n > 2 )
1072 {
1073 ny = TMath::Nint(TMath::Sqrt(n));
1074 nx = n/ny;
1075 }
1076
1077 TString stype(type);
1078 stype.ToUpper();
1079
1080 TString spectraName(Form("/%s/%s/%s/%s/%s-%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data(),flavour));
1081
1082 AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(spectraName.Data()));
1083
1084 AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
1085
1086 TObjArray* spectraBins(0x0);
1087 if ( spectra )
1088 {
1089 spectraBins = spectra->Bins();
1090 }
1091
1092 TCanvas* c = new TCanvas;
1093 c->Divide(nx,ny);
1094 c->Draw();
1095 gStyle->SetOptFit(1112);
1096
1097 c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
1098 (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
1099
1100
a58729a5 1101 TIter next(bins);
1102 AliAnalysisMuMuBinning::Range* r;
1afce1ce 1103 Int_t ci(0);
a58729a5 1104
1105 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1106 {
1107 TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
1108
1109 AliDebug(1,name.Data());
1afce1ce 1110
1111 AliAnalysisMuMuResult* spectraBin(0x0);
1112
1113 if ( spectraBins )
1114 {
1115 AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
1116
1117 spectraBin = sr->SubResult(subresultname);
1118
1119 AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
1120 }
1121
a58729a5 1122 TH1* h = MC()->Histo(name.Data());
1afce1ce 1123
1124 if ( spectraBin )
1125 {
1126 h = spectraBin->Minv();
1127 }
1128
a58729a5 1129 if (h)
1130 {
1afce1ce 1131 ++ci;
1132 c->cd(ci);
1133 gPad->SetLogy();
1134 if (xmin>0)
1135 {
1136 h->GetXaxis()->SetRangeUser(xmin,xmax);
1137 }
1138 h->Draw("histes");
1139
1140 TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
1141 if (f1)
1142 {
1143 f1->Draw("same");
1144 }
1145
1146 gPad->Modified();
1147 gPad->Update();
1148
1149 TObject* stats = h->FindObject("stats");
1150 if (stats)
1151 {
1152 stats->Draw("same");
1153 }
a58729a5 1154 }
1155 }
1156
1157 delete bins;
1158}
1159
1160//_____________________________________________________________________________
1afce1ce 1161void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
2754fd07 1162{
a58729a5 1163 /// Draw minv spectra for binning of given type
1164
1165 DrawMinv(type,particle,
1166 First(DimuonTriggerList()).Data(),
1167 First(EventSelectionList()).Data(),
1168 First(PairSelectionList()).Data(),
1afce1ce 1169 First(CentralitySelectionList()).Data(),
1170 subresultname,
1171 flavour);
1172}
1173
1174//___________________________________________________________________
1175void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
1176{
1177 // Actions in reponse to mouse button events.
1178
1179 TCanvas* c = static_cast<TCanvas*>(gTQSender);
1180 TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
1181 if (!pad) return;
1182
1183// if ((event == kButton1Down) ||
1184 if (event == kButton1Double)
1185 {
1186
1187// Float_t x = pad->AbsPixeltoX(px);
1188// Float_t y = pad->AbsPixeltoY(py);
1189// x = pad->PadtoX(x);
1190// y = pad->PadtoY(y);
1191
1192// std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
1193
1194 if ( sel && sel->InheritsFrom("TH1") )
1195 {
1196 TCanvas* clocal = new TCanvas;
1197 clocal->SetLogy();
1198 clocal->Draw();
1199 sel->Draw();
1200 }
1201 else
1202 {
1203 TList* list = pad->GetListOfPrimitives();
1204 TIter next(list);
1205 TObject* h;
1206
1207 while ( ( h = next() ) )
1208 {
1209 if ( h->InheritsFrom("TH1") )
1210 {
1211 TCanvas* clocal = new TCanvas;
1212 clocal->SetLogy();
1213 clocal->Draw();
1214 h->Draw();
1215 break;
1216 }
1217 }
1218
1219 }
1220
1221// std::cout << std::endl;
1222
1223 pad->Modified();
1224 }
1225
2754fd07 1226}
1227
2f331ac9 1228//_____________________________________________________________________________
1229TString
1230AliAnalysisMuMu::ExpandPathName(const char* file)
1231{
1232 // An expand method that lives alien URL as they are
1233 TString sfile;
1234
1235 if ( !sfile.BeginsWith("alien://") )
1236 {
1237 return gSystem->ExpandPathName(file);
1238 }
1239 else
1240 {
1241 if (!gGrid) TGrid::Connect("alien://");
1242 if (!gGrid) return "";
1243 }
1244
1245 return file;
1246}
1247
1248//_____________________________________________________________________________
1249TFile*
1250AliAnalysisMuMu::FileOpen(const char* file)
1251{
1252 // Open a file after expansion of its name
1253
1254 return TFile::Open(ExpandPathName(file).Data());
1255}
1256
a58729a5 1257//_____________________________________________________________________________
1258TString AliAnalysisMuMu::First(const TString& list) const
1259{
1260 TObjArray* a = list.Tokenize(",");
1261 if ( a->GetLast() < 0 ) return "";
1262
1263 TString rv = static_cast<TObjString*>(a->First())->String();
1264
1265 delete a;
1266
1267 return rv;
1268}
1269
1270//_____________________________________________________________________________
1271AliAnalysisMuMuSpectra*
1272AliAnalysisMuMu::FitParticle(const char* particle,
1273 const char* trigger,
1274 const char* eventType,
1275 const char* pairCut,
1276 const char* centrality,
1277 const AliAnalysisMuMuBinning& binning)
1278{
1279 // Fit the minv spectra to find the given particle
1280 // Returns an array of AliAnalysisMuMuResult objects
1281
1282 static int n(0);
1283
1284 TObjArray* bins = binning.CreateBinObjArray(particle);
1285 if (!bins)
1286 {
1287 AliError(Form("Did not get any bin for particle %s",particle));
1288 return 0x0;
1289 }
1290
1291 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1292 if ( !triggers->FindObject(trigger) )
1293 {
1294 AliDebug(1,Form("Did not find trigger %s",trigger));
1295 delete bins;
1296 delete triggers;
1297 return 0x0;
1298 }
1299 delete triggers;
1300
1301 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1302 if ( !events->FindObject(eventType) )
1303 {
1304 AliError(Form("Did not find eventType %s",eventType));
1305 delete bins;
1306 delete events;
1307 return 0x0;
1308 }
1309 delete events;
1310
1311 Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
1312
1313 if (ntrigger<=0)
1314 {
1315 AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
1316 delete bins;
1317 return 0x0;
1318 }
1319
1afce1ce 1320// binning.Print();
a58729a5 1321
1322 AliAnalysisMuMuSpectra* spectra(0x0);
1323
1324 AliAnalysisMuMuBinning::Range* bin;
1325 TIter next(bins);
1afce1ce 1326
1327 TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
1328 TIter nextFitType(fitTypeArray);
1329 TObjString* fitType;
1330 TString flavour;
a58729a5 1331
1332 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1333 {
1334 TString hname(Form("MinvUS%s",bin->AsString().Data()));
1335
a58729a5 1336 TH1* hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),hname.Data());
1337
1338 if (!hminv)
1339 {
1340 if (!fBinning && bin->IsNullObject() )
1341 {
1342 // old file, we only had MinvUSPt
1343 hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
1344 }
1345
1346 if (!hminv)
1347 {
1348 AliDebug(1,Form("Could not find histo %s",hname.Data()));
1349 continue;
1350 }
1351 }
1352
1353 hminv = static_cast<TH1*>(hminv->Clone(Form("minv%d",n++)));
1354
1355
1356 AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(*hminv,
1357 trigger,
1358 eventType,
1359 pairCut,
1360 centrality,
1361 *bin);
1362
1363 r->SetNofTriggers(ntrigger);
1364
1afce1ce 1365 nextFitType.Reset();
1366
1367 while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
a58729a5 1368 {
1afce1ce 1369 AliDebug(1,Form("<<<<<< fitType=%s bin=%s",fitType->String().Data(),bin->Flavour().Data()));
1370
1371 if ( fitType->String().BeginsWith("PSILOWMCTAILS") )
1372 {
1373 std::vector<Double_t> par;
1374 par = GetMCCB2Tails(*bin);
1375 if (!par.empty())
1376 {
1377 r->AddFit(fitType->String().Data(),par.size(),&par[0]);
1378 }
1379 }
1380 else
1381 {
1382 r->AddFit(fitType->String());
1383 }
a58729a5 1384 }
a58729a5 1385
1afce1ce 1386 flavour = bin->Flavour();
1387
a58729a5 1388 if (!spectra)
1389 {
1afce1ce 1390 TString spectraName(binning.GetName());
1391 if ( flavour.Length() > 0 )
1392 {
1393 spectraName += "-";
1394 spectraName += flavour;
1395 }
1396 spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
a58729a5 1397 }
1398
1399 spectra->AdoptResult(*bin,r);
1400
1401 if ( IsSimulation() )
1402 {
1403 SetNofInputParticles(*r);
1404 }
1405
1406
1407 } // loop on bins
1408
1afce1ce 1409 delete fitTypeArray;
a58729a5 1410 delete bins;
1411
1412 return spectra;
1413}
1414
1afce1ce 1415//_____________________________________________________________________________
1416std::vector<Double_t>
1417AliAnalysisMuMu::GetMCCB2Tails(const AliAnalysisMuMuBinning::Range& bin) const
1418{
1419 /// Get the tails from the associated simulation
1420
1421 std::vector<Double_t> par;
1422
1423 if (!SIM())
1424 {
1425 AliError("Cannot get MC tails without an associated simulation file !");
1426 return par;
1427 }
1428
1429
1430 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(SIM()->GetSpectra(bin.Type().Data(),bin.Flavour().Data()));
1431
1432 if (!s)
1433 {
1434 AliError(Form("Could not find spectra %s,%s for associated simulation",bin.Type().Data(),bin.Flavour().Data()));
1435 fAssociatedSimulation->MC()->Print("*:Ali*");
1436 return par;
1437 }
1438 else
1439 {
1440 AliDebug(1,Form("AliAnalysisMuMuSpectra* s = reinterpret_cast<AliAnalysisMuMuSpectra*>(%p)",s));
1441 }
1442
1443 AliAnalysisMuMuResult* r = s->GetResultForBin(bin);
1444
1445 if ( r )
1446 {
1447 AliAnalysisMuMuResult* r1 = r->SubResult("JPSI:1");
1448 if (r1)
1449 {
1450 TF1* func = static_cast<TF1*>(r1->Minv()->GetListOfFunctions()->FindObject("fitTotal"));
1451 if (func)
1452 {
1453 par.push_back(func->GetParameter("alphaLow"));
1454 par.push_back(func->GetParameter("nLow"));
1455 par.push_back(func->GetParameter("alphaUp"));
1456 par.push_back(func->GetParameter("nUp"));
1457 }
1458 }
1459 }
1460
1461 return par;
1462}
1463
2f331ac9 1464//_____________________________________________________________________________
1465ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
1466{
1467 // Get the expected (from OCDB scalers) trigger count
1468
1469 AliAnalysisTriggerScalers ts(runNumber,fgOCDBPath.Data());
1470
2754fd07 1471 TObjArray* triggers = TString(triggerList).Tokenize(",");
2f331ac9 1472 TObjString* trigger;
1473 TIter next(triggers);
1474 ULong64_t n(0);
1475
1476 while ( ( trigger = static_cast<TObjString*>(next()) ) )
1477 {
1478 AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
1479 if (item)
1480 {
1481 n += item->Value();
1482 }
1483 delete item;
1484 }
1485 delete triggers;
1486
1487 return n;
1488}
1489
1afce1ce 1490//_____________________________________________________________________________
1491AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
1492{
1493 /// get a given spectra
1494
1495 TString swhat(what);
1496 TString sflavour(flavour);
1497 swhat.ToUpper();
1498 sflavour.ToUpper();
1499
1500 TString spectraName(Form("/PSALL/%s/PP/%s/PSI-%s",
1501 First(fgDefaultDimuonTriggers).Data(),
1502 First(fgDefaultPairSelectionList).Data(),
1503 swhat.Data()));
1504
1505 if (sflavour.Length()>0)
1506 {
1507 spectraName += "-";
1508 spectraName += sflavour.Data();
1509 }
1510
1511 if (IsSimulation())
1512 {
1513 spectraName.ReplaceAll("PSALL",fgDefaultEventSelectionForSimulations.Data());
1514 spectraName.ReplaceAll(First(fgDefaultDimuonTriggers).Data(),fgDefaultDimuonTriggerForSimulations.Data());
1515 }
1516
1517 return SPECTRA(spectraName.Data());
1518}
1519
2f331ac9 1520//_____________________________________________________________________________
a58729a5 1521UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
1522 const char* eventSelection, Int_t runNumber)
2f331ac9 1523{
1524 TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
1525 TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
1526 TObjArray* a = TString(triggerList).Tokenize(" ");
1527 TIter next(a);
1528 TObjString* str;
1529
1530 UInt_t n(0);
1531
1532 TString sEventSelection(eventSelection);
1533 sEventSelection.ToUpper();
1534
1535 if ( kevent->FindObject(sEventSelection.Data()) )
1536 {
1537 while ( ( str = static_cast<TObjString*>(next()) ) )
1538 {
1539 if ( ktrigger->FindObject(str->String().Data()) )
1540 {
1541 if ( runNumber < 0 )
1542 {
1543 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
1544 }
1545 else
1546 {
1547 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
1548 }
1549 }
1550 }
1551 }
1552
1553 delete a;
1554 delete ktrigger;
1555 delete kevent;
1556 return n;
1557}
1558
1559//_____________________________________________________________________________
a58729a5 1560Bool_t
1561AliAnalysisMuMu::GetCollections(const char* rootfile,
1562 AliMergeableCollection*& mc,
1563 AliCounterCollection*& cc,
1564 AliAnalysisMuMuBinning*& bin,
1565 std::set<int>& runnumbers)
2f331ac9 1566{
a58729a5 1567 mc = 0x0;
1568 cc = 0x0;
1569 bin = 0x0;
2754fd07 1570
a58729a5 1571 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
1572 if (!f)
1573 {
1574 f = TFile::Open(rootfile);
1575 }
1576
1577 if ( !f || f->IsZombie() )
1578 {
1579 return kFALSE;
1580 }
1581
1582 f->GetObject("MC",mc);
1583 f->GetObject("CC",cc);
1584
1585 TIter next(f->GetListOfKeys());
1586 TKey* key;
1587
1588 while ( ( key = static_cast<TKey*>(next())) && !bin )
1589 {
1590 if ( strcmp(key->GetClassName(),"AliAnalysisMuMuBinning")==0 )
1591 {
1592 bin = dynamic_cast<AliAnalysisMuMuBinning*>(key->ReadObj());
1593 }
1594 }
2f331ac9 1595
a58729a5 1596 if (!mc || !cc)
2f331ac9 1597 {
a58729a5 1598 AliErrorClass("Old file. Please upgrade it!");
2f331ac9 1599
a58729a5 1600 return kFALSE;
2f331ac9 1601 }
a58729a5 1602
1603 // get run list
1604 TObjArray* runs = cc->GetKeyWords("run").Tokenize(",");
1605 runs->Sort();
1606 TIter nextRun(runs);
1607 TObjString* srun;
1608
1609 runnumbers.clear();
1610
1611 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
1612 {
1613 runnumbers.insert(srun->String().Atoi());
1614 }
1615
1616 delete runs;
1617
1618 return kTRUE;
1619}
1620
1621//_____________________________________________________________________________
1622Bool_t AliAnalysisMuMu::IsSimulation() const
1623{
1624 // whether or not we have MC information
1625
1afce1ce 1626 return ( fMergeableCollection->Histo(Form("/INPUT/%s/MinvUS",fgDefaultEventSelectionForSimulations.Data())) != 0x0 );
a58729a5 1627}
1628
1629//_____________________________________________________________________________
1630Int_t
1afce1ce 1631AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour)
a58729a5 1632{
1633 // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
1afce1ce 1634 // what="integrated" => fit only fully integrated MinvUS
a58729a5 1635 // what="pt" => fit MinvUS in pt bins
1636 // what="y" => fit MinvUS in y bins
1637 // what="pt,y" => fit MinvUS in (pt,y) bins
1638
1afce1ce 1639 TStopwatch timer;
1640
a58729a5 1641 if (!fMergeableCollection)
1642 {
1643 AliError("No mergeable collection. Consider Upgrade()");
1644 return 0;
1645 }
1646
1647 Int_t nfits(0);
1648
1649 TObjArray* triggerArray = fDimuonTriggers.Tokenize(",");
1650 TObjArray* eventTypeArray = fEventSelectionList.Tokenize(",");
1651 TObjArray* pairCutArray = fPairSelectionList.Tokenize(",");
1652 TObjArray* whatArray = TString(what).Tokenize(",");
1653
a58729a5 1654 TIter nextTrigger(triggerArray);
1655 TIter nextEventType(eventTypeArray);
1656 TIter nextPairCut(pairCutArray);
1657 TIter nextWhat(whatArray);
1658
1659 TObjString* trigger;
1660 TObjString* eventType;
1661 TObjString* pairCut;
1662 TObjString* swhat;
1663
1664 while ( ( swhat = static_cast<TObjString*>(nextWhat()) ) )
1665 {
1666 AliAnalysisMuMuBinning* binning(0x0);
1667
1668 if ( fBinning && swhat->String().Length() > 0 )
1669 {
1afce1ce 1670 binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
a58729a5 1671 }
1672 else
1673 {
1674 binning = new AliAnalysisMuMuBinning;
1675 binning->AddBin("psi",swhat->String().Data());
1676 }
2f331ac9 1677
a58729a5 1678 std::cout << "++++++++++++ swhat=" << swhat->String().Data() << std::endl;
2754fd07 1679
a58729a5 1680 if (!binning)
1681 {
1682 AliError("oups. binning is NULL");
1683 continue;
1684 }
1685 binning->Print();
1686
1687 nextTrigger.Reset();
2754fd07 1688
1689 while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
2f331ac9 1690 {
2754fd07 1691 AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
1692
1693 nextEventType.Reset();
1694
1695 while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
2f331ac9 1696 {
a58729a5 1697 AliDebug(1,Form("--EVENTTYPE %s",eventType->String().Data()));
2f331ac9 1698
2754fd07 1699 nextPairCut.Reset();
2f331ac9 1700
2754fd07 1701 while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
1702 {
a58729a5 1703 AliDebug(1,Form("----PAIRCUT %s",pairCut->String().Data()));
1704
1705 AliDebug(1,"----Fitting...");
1706
1707 AliAnalysisMuMuSpectra* spectra = FitParticle("psi",
1708 trigger->String().Data(),
1709 eventType->String().Data(),
1710 pairCut->String().Data(),
1711 "PP",
1712 *binning);
1713
1714 AliDebug(1,Form("----fitting done spectra = %p",spectra));
1715
1716 if ( spectra )
1717 {
1718 ++nfits;
1719
1720 TString id(Form("/%s/%s/PP/%s",eventType->String().Data(),
1721 trigger->String().Data(),
1722 pairCut->String().Data()));
1723
1724 TObject* o = fMergeableCollection->GetObject(id.Data(),spectra->GetName());
1725
1726 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
1727
1728 if (o)
1729 {
1730 AliWarning(Form("Replacing %s/%s",id.Data(),spectra->GetName()));
1731 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectra->GetName()));
1732 }
1733
1734 fMergeableCollection->Adopt(id.Data(),spectra);
1735
1736 spectra->Print();
1737 }
2754fd07 1738 }
2f331ac9 1739 }
1740 }
2f331ac9 1741 }
1742
a58729a5 1743 delete whatArray;
1744 delete triggerArray;
1745 delete eventTypeArray;
1746 delete pairCutArray;
1afce1ce 1747
1748 timer.Print();
1749
a58729a5 1750 if (nfits)
2f331ac9 1751 {
1afce1ce 1752 Update();
1753// ReOpen(fFilename,"UPDATE");
1754// fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
1755// ReOpen(fFilename,"READ");
a58729a5 1756 }
1afce1ce 1757
1758
a58729a5 1759 return nfits;
1760
1761}
1762
1763//_____________________________________________________________________________
1764void AliAnalysisMuMu::PlotBackgroundEvolution(const char* gfile, const char* triggerList, Double_t ymax, Bool_t fillBoundaries)
1765{
1766 // plot the graphs found in the file (and generated using the ComputeBackgroundEvolution() method)
1767
1768 TFile* f = TFile::Open(ExpandPathName(gfile).Data());
1769
1770 if ( !f || !f->IsOpen() )
1771 {
1772 return;
1773 }
1774
1775 SetColorScheme();
1776
1777
1778 TCanvas* c = new TCanvas("background-evolution","background-evolution");
1779
1780 c->Draw();
1781
1782 TLegend* l = new TLegend(0.4,0.6,0.97,0.97);
1783 l->SetFillColor(0);
1784 l->SetTextColor(AliAnalysisMuMu::kBlue);
1785 l->SetLineColor(AliAnalysisMuMu::kBlue);
1786
1787 TObjArray* triggers = TString(triggerList).Tokenize(",");
1788
1789 gStyle->SetOptTitle(0);
1790
1791 TObjString* str(0x0);
1792 TIter next(triggers);
1793 Int_t i(0);
1794 Int_t run1(99999999);
1795 Int_t run2(0);
1796
1797 std::set<int> runnumbers;
1798
1799 while ( ( str = static_cast<TObjString*>(next()) ) )
1800 {
1801 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
1802 if (!g) continue;
1803 for ( Int_t ir = 0; ir < g->GetN(); ++ir )
1804 {
1805 Int_t run = TMath::Nint(g->GetX()[ir]);
1806 runnumbers.insert(run);
1807 run1 = TMath::Min(run1,run);
1808 run2 = TMath::Max(run2,run);
1809 }
1810 }
1811
1812 AliInfoClass(Form("run1 %d run2 %d",run1,run2));
1813
1814 Double_t ymin(0);
1815
1816 TH2* hframe = new TH2F("hframe","hframe",run2-run1+1,run1,run2,100,ymin,ymax);
1817 hframe->Draw();
1818 hframe->GetXaxis()->SetNoExponent();
1819 hframe->GetYaxis()->SetTitle("Background percentage");
1820
1821 if (fillBoundaries)
1822 {
1823 AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath.Data());
1824 ts.DrawFills(ymin,ymax);
1825 }
1826
1827 next.Reset();
1828
1829 while ( ( str = static_cast<TObjString*>(next()) ) )
1830 {
1831 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
1832 if (!g)
1833 {
1834 AliErrorClass(Form("Graph mu_%s not found",str->String().Data()));
1835 continue;
1836 }
1837
1838 Int_t color(i+1);
1839
1840 if (i==0) color = AliAnalysisMuMu::kBlue;
1841 if (i==1) color = AliAnalysisMuMu::kOrange;
1842
1843 g->SetLineColor(color);
1844 g->SetMarkerColor(color);
1845 g->SetMarkerStyle(20+i);
1846
1847 g->Draw("LP");
1848
1849 TLegendEntry* le = l->AddEntry(g,str->String().Data(),"lp");
1850 le->SetTextColor(color);
1851
1852 g->GetYaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
1853 g->GetXaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
1854 // g->Print();
1855
1856 ++i;
2f331ac9 1857 }
2f331ac9 1858
a58729a5 1859 hframe->Draw("sameaxis");
1860
1861 l->Draw();
2f331ac9 1862}
1863
1864//_____________________________________________________________________________
2754fd07 1865void
1866AliAnalysisMuMu::PlotJpsiEvolution(const char* resultFile, const char* triggerList, Bool_t fillBoundaries,
a58729a5 1867 const char* efficiencyFile, Bool_t simulation)
2754fd07 1868{
a58729a5 1869 /// Will plot the Jpsi rate (or AccxEff if simulation=kTRUE) evolution.
1870 /// (JpsiRate is the number of Jpsi divided by the number of triggers)
1871
1872 std::map<int, float> efficiencies;
1873
2754fd07 1874 if ( efficiencyFile && strlen(efficiencyFile) > 0 )
1875 {
1876 std::ifstream in(gSystem->ExpandPathName(efficiencyFile));
1877 if (!in.bad())
1878 {
1879 char line[1024];
1880 int run;
a58729a5 1881 float eff;
1882 float dummy,errorl,errorh;
1883 int idummy;
2754fd07 1884 while ( in.getline(line,1023,'\n') )
1885 {
a58729a5 1886 sscanf(line,"%d, x[%d]=%f, y[%d]=%f, exl[%d]=%f, exh[%d]=%f, eyl[%d]=%f, eyh[%d]=%f",
1887 &run,&idummy,&dummy,&idummy,&eff,&idummy,&dummy,&idummy,&dummy,&idummy,&errorl,&idummy,&errorh);
1888
1889 AliDebugClass(1,Form("%09d %8.6f +- %8.6f ",run,eff,errorl+errorh));
1890
1891 efficiencies[run] = eff;
2754fd07 1892 }
1893 }
2754fd07 1894 }
2f331ac9 1895
2754fd07 1896 TFile* f = TFile::Open(gSystem->ExpandPathName(resultFile));
2f331ac9 1897
a58729a5 1898 std::set<int> runnumbers;
2f331ac9 1899
2754fd07 1900 TMap* m = static_cast<TMap*>(f->Get("rbr"));
1901
1902 TIter next(m);
2f331ac9 1903 TObjString* str;
2f331ac9 1904
2754fd07 1905 TObjArray files;
1906 files.SetOwner(kTRUE);
1907
1908 while ( ( str = static_cast<TObjString*>(next())) )
2f331ac9 1909 {
2754fd07 1910 files.Add(new TObjString(str->String()));
2f331ac9 1911 }
1912
2754fd07 1913 files.Sort();
2f331ac9 1914
2754fd07 1915 std::map<std::string, std::vector<float> > x_jpsirate;
1916 std::map<std::string, std::vector<float> > y_jpsirate;
1917 std::map<std::string, std::vector<float> > xerr_jpsirate;
1918 std::map<std::string, std::vector<float> > yerr_jpsirate;
2f331ac9 1919
2754fd07 1920 TIter nextTrigger(TString(triggerList).Tokenize(","));
1921 TObjString* trigger(0x0);
2f331ac9 1922
2754fd07 1923 int runMin(100000000);
1924 int runMax(0);
1925
1926 TIter nextFile(&files);
1927
a58729a5 1928 Double_t ymin(TMath::Limits<double>::Max());
1929 Double_t ymax(TMath::Limits<double>::Min());
1930
1931
2754fd07 1932 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2f331ac9 1933 {
a58729a5 1934 Double_t sumw(0);
1935 Double_t n(0);
1936
2754fd07 1937 TString triggerClass(trigger->String());
2f331ac9 1938
2754fd07 1939 nextFile.Reset();
a58729a5 1940
2754fd07 1941 while ( ( str = static_cast<TObjString*>(nextFile())) )
2f331ac9 1942 {
2754fd07 1943 TObjArray* a = static_cast<TObjArray*>(m->GetValue(str->String().Data()));
1944 if (!a) continue;
a58729a5 1945 AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(a->FindObject(triggerClass.Data()));
2754fd07 1946 if (!r) continue;
1947
1948 TString period;
1949 int aodtrain,esdpass,runnumber;
1950
1951 if ( DecodeFileName(str->String().Data(),period,esdpass,aodtrain,runnumber) )
2f331ac9 1952 {
a58729a5 1953 runnumbers.insert(runnumber);
1954
2754fd07 1955 runMin = TMath::Min(runMin,runnumber);
1956 runMax = TMath::Max(runMax,runnumber);
2f331ac9 1957
2754fd07 1958 x_jpsirate[triggerClass.Data()].push_back(runnumber);
1959 xerr_jpsirate[triggerClass.Data()].push_back(0.5);
2f331ac9 1960
2754fd07 1961 Double_t y(0.0);
1962 Double_t yerr(0.0);
a58729a5 1963 TString what("RateJpsi");
1964 if ( simulation )
1965 {
1966 what = "AccEffJpsi";
1967 }
2f331ac9 1968
a58729a5 1969 if ( TMath::Finite(r->GetValue("SigmaJpsi")) && r->NofTriggers() > 10 )
2f331ac9 1970 {
a58729a5 1971 y = 100*r->GetValue(what.Data());
1972 yerr = 100*r->GetErrorStat(what.Data());
2f331ac9 1973
a58729a5 1974 if (!efficiencies.empty() )
2754fd07 1975 {
a58729a5 1976 if (efficiencies.count(runnumber))
1977 {
1978 y /= ( efficiencies[runnumber] );
1979 }
1980 else
1981 {
1982 continue;
1983 }
2754fd07 1984 }
a58729a5 1985
1986 ymin = TMath::Min(ymin,y);
1987 ymax = TMath::Max(ymax,y);
1988
1989 sumw += y*r->NofTriggers();
1990 n += r->NofTriggers();
2f331ac9 1991 }
2754fd07 1992
1993 y_jpsirate[triggerClass.Data()].push_back(y);
1994 yerr_jpsirate[triggerClass.Data()].push_back(yerr);
2f331ac9 1995 }
2f331ac9 1996 }
a58729a5 1997
1998 AliInfoClass(Form("Trigger %30s ponderated mean is %7.2f",trigger->String().Data(),sumw/n));
2f331ac9 1999 }
2f331ac9 2000
2754fd07 2001 delete f;
2f331ac9 2002
a58729a5 2003 TString canvasName("cJpsiRateEvolution");
2004
2005 if ( !efficiencies.empty() )
2006 {
2007 canvasName += "Corr";
2008
2009 }
2010 TCanvas* c = new TCanvas(canvasName.Data(),canvasName.Data());
2f331ac9 2011
2754fd07 2012 c->Draw();
2f331ac9 2013
a58729a5 2014 Int_t nbins = runnumbers.size();
2015 Int_t xmin(runMin);
2016 Int_t xmax(runMax);
2017
2018 if ( CompactGraphs() )
2019 {
2020 xmin = 0;
2021 xmax = nbins-1;
2022 }
2f331ac9 2023
a58729a5 2024 TH2* h = new TH2F("h",Form("h;RunNumber;%s",(simulation ? "AccxEff (%)":"J/#psi per CMUL (%)")),
2025 nbins,xmin,xmax,100,ymin,ymax*1.2);
2f331ac9 2026
2754fd07 2027 gStyle->SetOptTitle(0);
2028 gStyle->SetOptStat(0);
2f331ac9 2029
a58729a5 2030 if (!CompactGraphs())
2031 {
2032 h->GetXaxis()->SetNoExponent();
2033 }
2034 else
2035 {
2036 std::set<int>::const_iterator it;
2037 int i(0);
2038
2039 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2040 {
2041 h->GetXaxis()->SetBinLabel(i,Form("%d",*it));
2042 ++i;
2043 }
2044 h->GetXaxis()->SetNdivisions(1,kFALSE);
2045
2046 }
2f331ac9 2047
2754fd07 2048 h->Draw();
2049
2050 if (fillBoundaries)
2051 {
a58729a5 2052 AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath);
2053 ts.DrawFills(ymin,ymax);
2754fd07 2054 }
2055
2056 h->Draw("sameaxis");
2f331ac9 2057
2754fd07 2058 //c->RedrawAxis("g");
2059
2060 nextTrigger.Reset();
2f331ac9 2061
2754fd07 2062 int i(0);
a58729a5 2063 int color[] = { 2,1,4,5,6 };
2754fd07 2064 int marker[] = { 20,23,25,21,22 };
2f331ac9 2065
2754fd07 2066 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2067 {
2068 std::vector<float>& x = x_jpsirate[trigger->String().Data()];
2069 std::vector<float>& y = y_jpsirate[trigger->String().Data()];
2070 std::vector<float>& xerr = xerr_jpsirate[trigger->String().Data()];
2071 std::vector<float>& yerr = yerr_jpsirate[trigger->String().Data()];
2072
2073 TGraphErrors* g = new TGraphErrors(x.size(),&x[0],&y[0],&xerr[0],&yerr[0]);
2074
a58729a5 2075 g->SetLineColor(1);//color[i]);
2754fd07 2076 g->SetMarkerColor(color[i]);
2077 g->SetMarkerStyle(marker[i]);
a58729a5 2078 g->SetMarkerSize(0.7);
2754fd07 2079 g->GetXaxis()->SetNoExponent();
a58729a5 2080
2081 if ( CompactGraphs() )
2082 {
2083 Compact(*g);
2084 }
2085
2086 g->Draw("P");
2087 TString gname(trigger->String());
2088 gname.ReplaceAll("-","_");
2089 g->SetName(gname.Data());
2754fd07 2090// g->Print();
2091
2092 Double_t m2 = g->GetMean(2);
2093
2094 TLine* line = new TLine(runMin,m2,runMax,m2);
2095 line->SetLineColor(color[i]);
2096 line->Draw();
2097
2098 AliInfoClass(Form("TRIGGER %s MEAN %7.2f",trigger->String().Data(),m2));
2099 ++i;
2100 }
2f331ac9 2101
2f331ac9 2102
2f331ac9 2103}
2104
2105//_____________________________________________________________________________
a58729a5 2106TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const char* event1,
2107 const char* trigger2, const char* event2,
2108 Bool_t drawFills,
2109 Bool_t asRejection) const
2f331ac9 2110{
a58729a5 2111 if (!CC()) return 0x0;
2112
2113 const std::set<int>& runnumbers = RunNumbers();
2114
2115 TGraphErrors* g = new TGraphErrors(runnumbers.size());
2116
2117 std::set<int>::const_iterator it;
2118 Int_t i(0);
2119
2120 Double_t ymin(TMath::Limits<double>::Max());
2121 Double_t ymax(TMath::Limits<double>::Min());
2122
2123 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2124 {
2125 Int_t runNumber = *it;
2126 Double_t n = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger1,event1,runNumber));
2127 Double_t d = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger2,event2,runNumber));
2128 if (n>0 && d>0)
2129 {
2130 Double_t y = n/d;
2131
2132 if ( fCorrectionPerRun )
2133 {
2134 Double_t xcorr,ycorr;
2135 fCorrectionPerRun->GetPoint(i,xcorr,ycorr); // note that the fact that xcorr==runNumber has been checked by the SetCorrectionPerRun method
1afce1ce 2136 y *= ycorr;
2137 // FIXME: should get the correction error here
a58729a5 2138 }
2139
2140 if ( asRejection ) y = 100*(1.0 - y);
2141 ymin = TMath::Min(ymin,y);
2142 ymax = TMath::Max(ymax,y);
2143 Double_t yerr = y*AliAnalysisMuMuResult::ErrorAB(n,TMath::Sqrt(n),d,TMath::Sqrt(d));
2144 g->SetPoint(i,runNumber,y);
2145 g->SetPointError(i,0.5,yerr);
2146
2147 ++i;
2148 }
2149
2150 }
2754fd07 2151
a58729a5 2152 TH2* hframe = new TH2F(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2153 Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2154 runnumbers.size()+50,
2155 *(runnumbers.begin())-25,
2156 *(runnumbers.rbegin())+25,100,0,ymax*1.3);
2f331ac9 2157
a58729a5 2158 gStyle->SetOptStat(0);
2159
2160 hframe->Draw();
2161
2162 hframe->GetXaxis()->SetNoExponent();
2163
2164 hframe->GetYaxis()->SetTitle(asRejection ? "Rejection (%)" : "Ratio");
2165
2166 g->Set(i);
2167 g->SetTitle(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2));
2168 g->GetXaxis()->SetNoExponent();
2169 g->Draw("lp");
2170
2171 AliAnalysisTriggerScalers ts(RunNumbers(),fgOCDBPath.Data());
2172
2173 if ( drawFills )
2754fd07 2174 {
a58729a5 2175 ts.DrawFills(ymin,ymax);
2176 g->Draw("lp");
2754fd07 2177 }
2f331ac9 2178
a58729a5 2179
2180 std::map<std::string, std::pair<int,int> > periods;
2181
2182 ts.GetLHCPeriodBoundaries(periods);
2183
2184 TLegend* legend = new TLegend(0.15,0.82,0.90,0.92);
2185 legend->SetFillColor(0);
2186 Int_t n(0);
2187
2188
2189 for ( std::map<std::string, std::pair<int,int> >::const_iterator pit = periods.begin(); pit != periods.end(); ++pit )
2190 {
2191 std::string period = pit->first;
2192 int run1 = (pit->second).first;
2193 int run2 = (pit->second).second;
2194 int nruns(0);
2195 for ( std::set<int>::const_iterator rit = RunNumbers().begin(); rit != RunNumbers().end(); ++ rit )
2196 {
2197 if ( (*rit) >= run1 && (*rit) <= run2 )
2198 {
2199 ++nruns;
2200 }
2201 }
2202 AliInfo(Form("Period %s runs %6d-%6d ; %d actual runs",period.c_str(),run1,run2,nruns));
2203
2204 g->Fit("pol0","+Q","",run1,run2);
2205 TF1* func = static_cast<TF1*>(g->GetListOfFunctions()->Last());
2206 if (func)
2207 {
2208 func->SetLineColor(2+n);
2209 legend->AddEntry(func,Form("%s %5.2f #pm %5.2f %s (rel. error %5.2f %%)",period.c_str(),func->GetParameter(0),func->GetParError(0),
2210 (asRejection ? "%":""),100*func->GetParError(0)/func->GetParameter(0)));
2211 ++n;
2212 }
2213 }
2214
2215 legend->SetNColumns(3);
2f331ac9 2216
a58729a5 2217 Double_t mean = TMath::Mean(g->GetN(),g->GetY());
2218 Double_t rms = TMath::RMS(g->GetN(),g->GetY());
2f331ac9 2219
a58729a5 2220 legend->AddEntry("",Form("Mean %5.2f RMS %5.2f (%5.2f %%)",mean,rms,(mean) ? 100.0*rms/mean : 0.0),"");
2f331ac9 2221
a58729a5 2222 legend->Draw();
2223
2224 return g;
2f331ac9 2225}
2226
2227//_____________________________________________________________________________
a58729a5 2228void AliAnalysisMuMu::Print(Option_t* opt) const
2f331ac9 2229{
a58729a5 2230 /// printout
2231 std::cout << "Reading from file : " << fFilename.Data() << std::endl;
2232 std::cout << "List of dimuon triggers to consider : " << DimuonTriggerList() << std::endl;
2233 std::cout << "List of muon triggers to consider : " << MuonTriggerList() << std::endl;
2234 std::cout << "List of MB triggers to consider : " << MinbiasTriggerList() << std::endl;
2235 std::cout << "Event selection list : " << EventSelectionList() << std::endl;
2236 std::cout << "Pair selection list : " << PairSelectionList() << std::endl;
2f331ac9 2237
a58729a5 2238 std::cout << RunNumbers().size() << " runs";
2239 if ( fCorrectionPerRun )
2f331ac9 2240 {
a58729a5 2241 std::cout << " with correction factors";
2f331ac9 2242 }
a58729a5 2243 std::cout << std::endl;
2244 Int_t i(0);
2245 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2f331ac9 2246 {
a58729a5 2247 std::cout << (*it);
2248 if ( fCorrectionPerRun )
2249 {
2250 std::cout << Form("(%e)",fCorrectionPerRun->GetY()[i]);
2251 }
2252 std::cout << ",";
2253 ++i;
2f331ac9 2254 }
a58729a5 2255 std::cout << std::endl;
2f331ac9 2256
a58729a5 2257 TString sopt(opt);
2258 sopt.ToUpper();
2f331ac9 2259
a58729a5 2260 if ( sopt.Contains("BIN") && BIN() )
2261 {
2262 std::cout << "Binning : " << std::endl;
2263 TString topt(sopt);
2264 topt.ReplaceAll("BIN","");
2265 BIN()->Print(topt.Data());
2266 }
2267 if ( sopt.Contains("MC") && MC() )
2268 {
2269 TString topt(sopt);
2270 topt.ReplaceAll("MC","");
2271 MC()->Print(topt.Data());
2272 }
2273 if ( sopt.Contains("CC") && CC() )
2274 {
2275 CC()->Print("trigger/event");
2276 }
2f331ac9 2277
a58729a5 2278 if ( sopt.Contains("SIZE") )
2279 {
2280 TFile* f = ReOpen(fFilename,"READ");
2281 TIter next(f->GetListOfKeys());
2282 TKey* key;
2283
2284 while ( ( key = static_cast<TKey*>(next()) ) )
2285 {
2286 std::cout << key->GetName() << " " << key->GetNbytes() << " " << key->GetObjlen() << std::endl;
2287 }
2288 }
2f331ac9 2289}
2290
2291//_____________________________________________________________________________
a58729a5 2292TObjArray*
2293AliAnalysisMuMu::ReadFileList(const char* filelist)
2f331ac9 2294{
a58729a5 2295 //
2296 // read the filelist and try to order it by runnumber
2297 //
2298 // filelist can either be a real filelist (i.e. a text file containing
2299 // root filenames) or a root file itself.
2300 //
2301
2302 char line[1024];
2303
2304 TObjArray* files = new TObjArray;
2305 files->SetOwner(kTRUE);
2f331ac9 2306
a58729a5 2307 TString sfilelist(ExpandPathName(filelist));
2f331ac9 2308
a58729a5 2309 if ( sfilelist.EndsWith(".root") )
2f331ac9 2310 {
a58729a5 2311 files->Add(new TObjString(sfilelist.Data()));
2312 return files;
2f331ac9 2313 }
2314
a58729a5 2315 std::set<int> runnumbers;
2316 std::map<int,std::string> filemap;
2f331ac9 2317
a58729a5 2318 std::ifstream in(sfilelist.Data());
2f331ac9 2319
a58729a5 2320 TString period;
2321 int aodtrain,esdpass,runnumber;
2f331ac9 2322
a58729a5 2323 while ( in.getline(line,1022,'\n') )
2324 {
2325 DecodeFileName(line,period,esdpass,aodtrain,runnumber);
2326
2327 AliDebugClass(1,Form("line %s => period %s esdpass %d aodtrain %d runnumber %09d",
2328 line,period.Data(),esdpass,aodtrain,runnumber));
2329
2330 filemap.insert(std::make_pair<int,std::string>(runnumber,line));
2331 runnumbers.insert(runnumber);
2332 }
2f331ac9 2333
a58729a5 2334 in.close();
2f331ac9 2335
a58729a5 2336 std::set<int>::const_iterator it;
2f331ac9 2337
a58729a5 2338 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2f331ac9 2339 {
a58729a5 2340 files->Add(new TObjString(filemap[*it].c_str()));
2f331ac9 2341 }
2342
a58729a5 2343 return files;
2344}
2f331ac9 2345
a58729a5 2346//_____________________________________________________________________________
2347TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
2348{
2349 /// Tries to reopen the file with a new mode
2f331ac9 2350
a58729a5 2351 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(filename));
2352
2353 if (f)
2f331ac9 2354 {
a58729a5 2355 delete f;
2f331ac9 2356 }
2357
1afce1ce 2358 f = TFile::Open(filename,mode);
a58729a5 2359
2360 if ( !f || !f->IsOpen() )
2361 {
2362 AliError(Form("Cannot open file %s in mode %s",filename,mode));
2363 return 0x0;
2364 }
2365
2366 return f;
2367}
2368
2369//_____________________________________________________________________________
2370Int_t
2371AliAnalysisMuMu::RunNumberFromFileName(const char* filename)
2372{
2373 TString period;
2374 int esdpass,aodtrain,runnumber;
2375 Bool_t ok = DecodeFileName(filename,period,esdpass,aodtrain,runnumber);
2376 if ( ok ) return runnumber;
2377 return -1;
2f331ac9 2378}
2379
2380//_____________________________________________________________________________
2381void AliAnalysisMuMu::SetColorScheme()
2382{
2383 new TColor(AliAnalysisMuMu::kBlue,4/255.0,44/255.0,87/255.0,"my blue");
2384 new TColor(AliAnalysisMuMu::kOrange,255/255.0,83/255.0,8/255.0,"my orange");
2385 new TColor(AliAnalysisMuMu::kGreen,152/255.0,202/255.0,52/255.0,"my green");
2386
2387 gStyle->SetGridColor(AliAnalysisMuMu::kBlue);
2388
2389 gStyle->SetFrameLineColor(AliAnalysisMuMu::kBlue);
2390 gStyle->SetAxisColor(AliAnalysisMuMu::kBlue,"xyz");
2391 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue,"xyz");
2392
2393 gStyle->SetTitleColor(AliAnalysisMuMu::kBlue);
2394 gStyle->SetTitleTextColor(AliAnalysisMuMu::kBlue);
2395 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue);
2396 gStyle->SetStatTextColor(AliAnalysisMuMu::kBlue);
2397
2398 gStyle->SetOptStat(0);
2399}
2400
2401//_____________________________________________________________________________
1afce1ce 2402Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
2f331ac9 2403{
a58729a5 2404 /// Sets the graph used to correct values per run
2405 delete fCorrectionPerRun;
2406 fCorrectionPerRun=0x0;
2f331ac9 2407
a58729a5 2408 // check that corr has the same runs as we do
2f331ac9 2409
a58729a5 2410 Int_t i(0);
2f331ac9 2411
a58729a5 2412 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2f331ac9 2413 {
a58729a5 2414 Int_t corrRun = TMath::Nint(corr.GetX()[i]);
2415 if (corrRun != *it)
2f331ac9 2416 {
a58729a5 2417 AliError(Form("%d-th run mistmatch %d vs %d",i,corrRun,*it));
2418
2419 return kFALSE;
2f331ac9 2420 }
a58729a5 2421 ++i;
2f331ac9 2422 }
2423
1afce1ce 2424 fCorrectionPerRun = new TGraphErrors(corr.GetN());
2425
2426 TFormula* tformula(0x0);
2427 if ( strlen(formula) > 0 )
2428 {
2429 tformula = new TFormula("SetCorrectionPerRunFormula",formula);
2430 }
2431
2432 i = 0;
2433
2434 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2435 {
2436 Double_t y = corr.GetY()[i];
2437
2438 if ( tformula )
2439 {
2440 y = tformula->Eval(y);
2441 }
2442 fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
2443 ++i;
2444 }
a58729a5 2445
1afce1ce 2446 delete formula;
2447
a58729a5 2448 return kTRUE;
2449}
2450
2451//_____________________________________________________________________________
2452void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuResult& r)
2453{
2454 /// Set the "NofInput" variable(s) of one result
2f331ac9 2455
a58729a5 2456 TString hname(Form("MinvUS%s",r.Bin().AsString().Data()));
2457
2458 TH1* hinput = fMergeableCollection->Histo("/INPUT/INYRANGE",hname.Data());
2459
2460 if (!hinput)
2461 {
2462 AliError(Form("Got a simulation file where I did not find histogram /INPUT/INYRANGE/%s",hname.Data()));
1afce1ce 2463
a58729a5 2464 }
2465 else
2466 {
2467 r.SetNofInputParticles(*hinput);
2468 }
2f331ac9 2469}
2470
1afce1ce 2471//_____________________________________________________________________________
2472AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
2473{
2474 /// Shortcut method to get to a spectra
2475 if (!MC()) return 0x0;
2476
2477 return static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(fullpath));
2478}
2479
2f331ac9 2480//_____________________________________________________________________________
a58729a5 2481void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
2482 Bool_t compact,
2483 Bool_t orderByTriggerCount)
2f331ac9 2484{
2485 // Give the fraction of triggers (in triggerList) relative
2486 // to what is expected in the scalers
2487
2488 TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
2489
2490 AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
2491
2492 AliLog::SetGlobalLogLevel(AliLog::kFatal);
2493
a58729a5 2494 if (!fMergeableCollection || !fCounterCollection) return;
2f331ac9 2495
2496 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2497 TIter nextRun(runs);
2498
2499 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2500 TIter nextTrigger(triggers);
2501
2502 TObjString* srun;
2503 TObjString* strigger;
2504
2505 TString striggerList(triggerList);
2506
2507 ULong64_t total(0);
2508 ULong64_t totalExpected(0);
a58729a5 2509 TString msg;
2510 std::multimap<ULong64_t,std::string> messages;
2f331ac9 2511
2512 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2513 {
a58729a5 2514 msg.Form("RUN %09d ",srun->String().Atoi());
2754fd07 2515
a58729a5 2516 if (!compact)
2517 {
2518 msg += "\n";
2519 }
2520
2521 ULong64_t nmax(0);
2522
2f331ac9 2523 nextTrigger.Reset();
2524
2525 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2526 {
2527 if ( !striggerList.Contains(strigger->String().Data()) )
2528 {
2529 continue;
2530 }
a58729a5 2531 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2532 strigger->String().Data(),"ALL",srun->String().Atoi())));
2f331ac9 2533
2534 ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
2535
2536
a58729a5 2537 nmax = TMath::Max(n,nmax);
2538
2f331ac9 2539 total += n;
2540 totalExpected += expected;
2541
1afce1ce 2542 msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
2543 (n>expected ? "!" : " "));
2f331ac9 2544
2545 if ( expected > 0 ) {
a58729a5 2546 msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
2f331ac9 2547 }
1afce1ce 2548
2754fd07 2549 if (!compact)
2550 {
a58729a5 2551 msg += "\n";
2552 }
2553 }
2554 if (nmax>0)
2555 {
2556 if (!orderByTriggerCount)
2557 {
2558 std::cout << msg.Data() << std::endl;
2559 }
2560 else
2561 {
2562 messages.insert(std::make_pair<ULong64_t,std::string>(nmax,msg.Data()));
2754fd07 2563 }
2f331ac9 2564 }
2f331ac9 2565 }
2566
a58729a5 2567 std::multimap<ULong64_t,std::string>::const_reverse_iterator it;
2568
2569 ULong64_t current(0);
2570 Int_t n(0);
2571
2572 for ( it = messages.rbegin(); it != messages.rend(); ++it )
2573 {
2574 ++n;
2575 current += it->first;
1afce1ce 2576 Double_t percent = ( total > 0.0 ? current*100.0/total : 0.0);
2577 std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",percent,n ) << std::endl;
a58729a5 2578 }
2579
2580 std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
2581 total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
2f331ac9 2582
a58729a5 2583
2f331ac9 2584 AliLog::SetGlobalLogLevel(oldLevel);
09d5920f 2585 delete triggers;
2586 delete runs;
2f331ac9 2587}
2588
2589//_____________________________________________________________________________
a58729a5 2590void AliAnalysisMuMu::UnsetCorrectionPerRun()
2591{
2592 // drop the correction factors
2593 delete fCorrectionPerRun;
2594 fCorrectionPerRun=0x0;
2595}
2596
1afce1ce 2597//_____________________________________________________________________________
2598void AliAnalysisMuMu::Update()
2599{
2600 /// update the current file with memory
2601
2602 ReOpen(fFilename,"UPDATE");
2603
2604 MC()->Write("MC",TObject::kSingleKey);
2605
2606 ReOpen(fFilename,"READ");
2607}
2608
a58729a5 2609//_____________________________________________________________________________
2610Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
2611{
2612 /// Upgrade a file
2613 AliAnalysisMuMu m(filename);
2614
2615 return m.Upgrade();
2616}
2617
2618//_____________________________________________________________________________
2619Bool_t AliAnalysisMuMu::Upgrade()
2620{
2621 /// Upgrade the current file
2622 /// - from single list to one key per object, if needed
2623 /// - from histogramCollection to mergeableCollection, if needed
2624
2625 TFile* f = ReOpen(fFilename,"UPDATE");
2626
2627 TList* list = static_cast<TList*>(f->Get("chist"));
2628
2629 if (list)
2630 {
2631 // really old file where everything was in a single list
2632
2633 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(list->At(0));
2634 AliCounterCollection* cc = static_cast<AliCounterCollection*>(list->At(1));
2635
2636 AliMergeableCollection* mc = hc->Convert();
2637
2638 f->cd();
2639
2640 mc->Write("MC",TObject::kSingleKey);
2641 cc->Write("CC",TObject::kSingleKey);
2642
2643 f->Delete("chist;*");
2644
2645 f->Write();
2646
2647 }
2648 else
2649 {
2650 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(f->Get("HC"));
2651
2652 if ( hc )
2653 {
2654 // old file with histogram collection instead of mergeable collection
2655
2656 AliMergeableCollection* mc = hc->Convert();
2657
2658 f->cd();
2659
2660 mc->Write("MC",TObject::kSingleKey);
2661
2662 f->Delete("HC;*");
2663
2664 f->Write();
2665 }
2666 }
2667
2668 delete f;
2669
2670 return kTRUE;
2671}
2672
2673//_____________________________________________________________________________
1afce1ce 2674AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
2f331ac9 2675{
1afce1ce 2676 /// Correct one spectra
a58729a5 2677
1afce1ce 2678 if (!SIM())
2679 {
2680 AliError("Cannot compute corrected yield without associated MC file !");
2681 return 0x0;
2682 }
2683
2684 const char* accEffSubResultName="COUNTJPSI:1";
a58729a5 2685
1afce1ce 2686 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
2687 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
a58729a5 2688
1afce1ce 2689 if ( !realSpectra )
2690 {
2691 AliError("could not get real spectra");
2692 return 0x0;
2693 }
a58729a5 2694
1afce1ce 2695 if ( !simSpectra)
a58729a5 2696 {
1afce1ce 2697 AliError("could not get sim spectra");
a58729a5 2698 return 0x0;
2699 }
2700
1afce1ce 2701 realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
2702
2703 Update();
2704
2705 return realSpectra;
2706}
2707
2708//_____________________________________________________________________________
2709AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
2710{
2711 if (!SIM())
2712 {
2713 AliError("Cannot compute corrected yield without associated MC file !");
2714 return 0x0;
2715 }
2716
2717 const char* accEffSubResultName="COUNTJPSI:1";
2718
2719 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
2720
2721 if ( !realSpectra )
2722 {
2723 AliError("could not get real spectra");
2724 return 0x0;
2725 }
2726
2727 if (!realSpectra->HasValue("CoffNofJpsi"))
2728 {
2729 if (!CorrectSpectra(type,flavour))
2730 {
2731 AliError("Could not get corrected spectra");
2732 return 0x0;
2733 }
2734 }
2735
2736 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
2737
2738 if ( !simSpectra)
2f331ac9 2739 {
a58729a5 2740 AliErrorClass("could not get sim spectra");
2741 return 0x0;
2f331ac9 2742 }
a58729a5 2743
1afce1ce 2744 Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
2745 Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
2746 Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
a58729a5 2747
1afce1ce 2748 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
2749 TObjArray* bins = binning->CreateBinObjArray();
2750 TIter nextBin(bins);
2751 AliAnalysisMuMuBinning::Range* bin;
2752 Int_t i(0);
2753 AliAnalysisMuMuResult* r;
a58729a5 2754
1afce1ce 2755 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
2756 {
2757 r = static_cast<AliAnalysisMuMuResult*>(realSpectra->Bins()->At(i));
2758
2759 StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
2760
2761 AliAnalysisMuMuResult* rsim = static_cast<AliAnalysisMuMuResult*>(simSpectra->Bins()->At(i));
2762
2763 Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
2764 Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2765 nofCINT7,TMath::Sqrt(nofCINT7),
2766 nofCMUL7,TMath::Sqrt(nofCMUL7));
2767
2768 r->Set("MBR",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2769 nofCINT7,TMath::Sqrt(nofCINT7)));
2770
2771 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
2772
2773 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
2774 mbeq,mbeqError);
2775
2776 r->Set("YJpsi",yield,yieldError);
2777
2778 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
2779 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
2780
2781 ++i;
2782 }
a58729a5 2783
1afce1ce 2784 delete bins;
2785
2786 Update();
2787
2788 return realSpectra;
2f331ac9 2789}
2790
1afce1ce 2791////_____________________________________________________________________________
2792//AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* realFile, const char* simFile,
2793// const char* type)
2794//{
2795// const char* accEffSubResultName="COUNTJPSI-1";
2796//
2797// AliAnalysisMuMu real(realFile);
2798// AliAnalysisMuMu sim(simFile);
2799//
2800//
2801// AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(real.MC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
2802// AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(sim.MC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
2803//
2804// if ( !realSpectra )
2805// {
2806// AliErrorClass("could not get real spectra");
2807// return 0x0;
2808// }
2809//
2810// if ( !simSpectra)
2811// {
2812// AliErrorClass("could not get sim spectra");
2813// return 0x0;
2814// }
2815//
2816// AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
2817// corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
2818//
2819// Double_t nofCMUL7 = real.CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
2820// Double_t nofCINT7 = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
2821// Double_t nofCINT7w0MUL = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
2822//
2823// AliAnalysisMuMuBinning* binning = corrSpectra->Binning();
2824// TObjArray* bins = binning->CreateBinObjArray();
2825// TIter nextBin(bins);
2826// AliAnalysisMuMuBinning::Range* bin;
2827// Int_t i(0);
2828// AliAnalysisMuMuResult* r;
2829//
2830// while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
2831// {
2832// r = static_cast<AliAnalysisMuMuResult*>(corrSpectra->Bins()->At(i));
2833//
2834// AliAnalysisMuMuResult* rsim = static_cast<AliAnalysisMuMuResult*>(simSpectra->Bins()->At(i));
2835//
2836// Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
2837// Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2838// nofCINT7,TMath::Sqrt(nofCINT7),
2839// nofCMUL7,TMath::Sqrt(nofCMUL7));
2840//
2841// r->Set("MBR",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2842// nofCINT7,TMath::Sqrt(nofCINT7)));
2843//
2844// Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
2845//
2846// Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
2847// mbeq,mbeqError);
2848//
2849// r->Set("YJpsi",yield,yieldError);
2850//
2851// r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
2852// r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
2853//
2854// ++i;
2855// }
2856//
2857// delete bins;
2858//
2859// return corrSpectra;
2860//}
2861
2f331ac9 2862//_____________________________________________________________________________
a58729a5 2863AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* realFile, const char* simFile, const char* type,
2864 const char* direction)
2f331ac9 2865{
a58729a5 2866 /// Compute the RAB...
2867 Double_t rapidityShift = 0.465;// 0.5*TMath::Log(208.0/82.0);
2868 const Double_t sqrts=5.023;
2869 const Double_t ymax=TMath::Log(sqrts*1000.0/3.096916);
2870 const Double_t tab = 0.093e-6; // nb^-1
2871 const Double_t tabError = 0.0035E-6; // nb^-1
1afce1ce 2872 const char* accEffSubResultName="COUNTJPSI:1";
2f331ac9 2873
a58729a5 2874 TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
2875 ydist.SetParameter(0,1.);
2876
2877 //Normalization to the values presented by Zaida and Rosana on January 11th 2013 https://indico.cern.ch/conferenceDisplay.py?confId=224985 slide 22
2878 // Normalization is done in the rapidity range 2.75<y<3.25 where Rosanas values is 230.8+212.1
2879 Double_t y1_norma= 2.75/ymax;
2880 Double_t y2_norma= 3.25/ymax;
2881 Double_t normalization = 0.25*(230.8+212.1)/ydist.Integral(y1_norma, y2_norma);
2882 ydist.SetParameter(0,normalization);
2883// AliInfoClass(Form("ymax=%e normalization=%f",ymax,ydist.Integral(y1_norma, y2_norma)));
2f331ac9 2884
a58729a5 2885 AliAnalysisMuMu real(realFile);
2886 AliAnalysisMuMu sim(simFile);
2f331ac9 2887
2f331ac9 2888
a58729a5 2889 AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(real.MC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
2890 AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(sim.MC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
2f331ac9 2891
a58729a5 2892 if ( !realSpectra )
2893 {
2894 AliErrorClass("could not get real spectra");
2895 return 0x0;
2896 }
2f331ac9 2897
a58729a5 2898 if ( !simSpectra)
2899 {
2900 AliErrorClass("could not get sim spectra");
2901 return 0x0;
2902 }
2903
2904 AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
2905 corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
2906
2907 Double_t nofCMUL7 = real.CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
2908 Double_t nofCINT7 = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
2909 Double_t nofCINT7w0MUL = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
2910
2911 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
2912 TObjArray* bins = binning->CreateBinObjArray();
2913 TIter nextBin(bins);
2914 AliAnalysisMuMuBinning::Range* bin;
2915 Int_t i(0);
2916 AliAnalysisMuMuResult* r;
2917
2918 Int_t n = bins->GetLast();
2919
2920 TObjArray finalBins(n+1);
2921 finalBins.SetOwner(kTRUE);
2f331ac9 2922
a58729a5 2923 TObjArray finalResults(n+1);
2924 finalResults.SetOwner(kFALSE);
2925
2926 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
2f331ac9 2927 {
a58729a5 2928 Double_t ylowlab = bin->Xmin();
2929 Double_t yhighlab = bin->Xmax();
2930
2931 Double_t ylowcms, yhighcms;
2932 Double_t ylownorm, yhighnorm;
2f331ac9 2933
a58729a5 2934 if ( bin->IsNullObject() )
2935 {
2936 ylowlab = -4;
2937 yhighlab = -2.5;
2938 }
2f331ac9 2939
a58729a5 2940 if ( strcmp(direction,"pPb")==0 )
2f331ac9 2941 {
a58729a5 2942 ylowcms = TMath::Abs(yhighlab) - rapidityShift;
2943 yhighcms = TMath::Abs(ylowlab) - rapidityShift;
2944 ylownorm = ylowcms/ymax;
2945 yhighnorm = yhighcms/ymax;
2f331ac9 2946 }
2947 else
2948 {
a58729a5 2949 ylowcms = ylowlab - rapidityShift;
2950 yhighcms = yhighlab - rapidityShift;
2951 ylownorm = -yhighcms/ymax;
2952 yhighnorm = -ylowcms/ymax;
2f331ac9 2953 }
2954
2f331ac9 2955
a58729a5 2956 Double_t brsigmapp = ydist.Integral(ylownorm,yhighnorm);
2957 Double_t brsigmappError = 0.0; // FIXME
2958
2959 AliInfoClass(Form("y range : LAB %f ; %f CMS %f ; %f -> ynorm : %f ; %f -> BR x sigmapp = %f",
2960 ylowlab,yhighlab,ylowcms,yhighcms,ylownorm,yhighnorm,brsigmapp));
2961
2962 r = static_cast<AliAnalysisMuMuResult*>(corrSpectra->Bins()->At(i)->Clone());
2963
2964 AliAnalysisMuMuResult* rsim = static_cast<AliAnalysisMuMuResult*>(simSpectra->Bins()->At(i));
2965
2966 Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
2967 Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2968 nofCINT7,TMath::Sqrt(nofCINT7),
2969 nofCMUL7,TMath::Sqrt(nofCMUL7));
2970
2971 r->Set("MBR",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
2972 nofCINT7,TMath::Sqrt(nofCINT7)));
2973
2974 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
2975
2976 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
2977 mbeq,mbeqError);
2978
2979 r->Set(Form("Y%sJpsi",direction),yield,yieldError);
2980
2981 Double_t raa = yield/(tab*brsigmapp);
2982 Double_t raaError = AliAnalysisMuMuResult::ErrorABC(yield,yieldError,
2983 tab,tabError,
2984 brsigmapp,brsigmappError);
2985 r->Set(Form("R%sJpsi",direction),raa,raaError);
2986
2987 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
2988 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
2989
2990 AliAnalysisMuMuBinning::Range* bincm = new AliAnalysisMuMuBinning::Range(bin->Particle(),bin->Type(),ylowcms,yhighcms);
2991
2992 r->SetBin(*bincm);
2993
2994 finalBins.Add(bincm);
2995 finalResults.Add(r);
2996
2997 ++i;
2f331ac9 2998 }
2999
a58729a5 3000 delete bins;
2f331ac9 3001
a58729a5 3002 AliAnalysisMuMuSpectra* spectra = new AliAnalysisMuMuSpectra(type,direction);
3003
3004 for ( i = 0; i <= n; ++i )
3005 {
3006 Int_t j(i);
3007 if ( strcmp(direction,"pPb")==0 )
3008 {
3009 j = n-i;
3010 }
3011
3012 r = static_cast<AliAnalysisMuMuResult*>(finalResults.At(j));
3013
3014 bin = static_cast<AliAnalysisMuMuBinning::Range*>(finalBins.At(j));
3015
3016 spectra->AdoptResult(*bin,r);
3017 }
3018
3019
3020 delete corrSpectra;
2f331ac9 3021
a58729a5 3022 return spectra;
2f331ac9 3023}
3024
3025//_____________________________________________________________________________
a58729a5 3026TGraph* AliAnalysisMuMu::ResultEvolution(const char* runlist, const char* period, const char* what, Bool_t forceRecomputation)
2f331ac9 3027{
a58729a5 3028 std::vector<int> runs;
3029 AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3030
3031 TGraphErrors* g = new TGraphErrors(runs.size());
2f331ac9 3032
a58729a5 3033 TString direction("Pbp");
2f331ac9 3034
a58729a5 3035 if (TString(period) == "LHC13b" ||
3036 TString(period) == "LHC13c" ||
3037 TString(period) == "LHC13d" ||
3038 TString(period) == "LHC13e"
3039 )
3040 {
3041 direction = "pPb";
3042 }
2f331ac9 3043
a58729a5 3044 AliInfoClass(Form("period %s direction %s",period,direction.Data()));
2f331ac9 3045
a58729a5 3046 Double_t weightedMean(0.0);
3047 Double_t sumOfWeights(0.0);
3048
3049 Double_t mean(0.0);
3050 TString subResultName("");
3051 TString swhat(what);
3052
3053 for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
2f331ac9 3054 {
a58729a5 3055 Int_t runNumber = runs[i];
3056
3057 AliInfoClass(Form("RUN %09d",runNumber));
3058
3059 TString realFile(Form("RUNBYRUN/%s_muon_calo_AODMUON000_%09d.saf.root",period,runNumber));
3060
3061 TString simFileName(Form("RUNBYRUN/SIM_JPSI_%s_pp503_newalign_%09d.saf.root",period,runNumber));
3062 if ( direction == "pPb" )
3063 {
3064 simFileName = Form("RUNBYRUN/SIM_JPSI_%s_pp503_%09d.saf.root",period,runNumber);
3065 }
3066
3067 TString simFile(simFileName);
3068
3069 TString resultName(Form("%s%sJpsi",what,direction.Data()));
3070
3071 if ( swhat == "MBR")
3072 {
3073 resultName = "MBR";
3074 }
3075 else if ( swhat.Contains("Acc") )
3076 {
3077 resultName.ReplaceAll(direction,"");
3078 }
3079
3080 AliAnalysisMuMu mreal(realFile);
3081
3082 AliAnalysisMuMuSpectra* real = static_cast<AliAnalysisMuMuSpectra*>(mreal.MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-"));
3083
3084 if (!real || forceRecomputation)
3085 {
3086 mreal.Jpsi();
3087 real = static_cast<AliAnalysisMuMuSpectra*>(mreal.MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-"));
3088 if (!real)
3089 {
3090 AliErrorClass(Form("Could not get real spectra for run %d",runNumber));
3091 return 0x0;
3092 }
3093 }
3094
3095 AliAnalysisMuMu msim(simFile);
3096
3097 AliAnalysisMuMuSpectra* sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-"));
3098
3099 if (!sim || forceRecomputation)
3100 {
3101 msim.SetEventSelectionList("ALL");
3102 msim.Jpsi();
3103 sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-"));
3104 if (!sim)
3105 {
3106 AliErrorClass(Form("Could not get sim spectra for run %d",runNumber));
3107 return 0x0;
3108 }
3109 }
3110
3111 AliAnalysisMuMuSpectra* corrected = AliAnalysisMuMu::RABy(realFile.Data(),simFile.Data(),"",direction.Data());
3112
3113 AliAnalysisMuMuResult* result = static_cast<AliAnalysisMuMuResult*>(corrected->Bins()->First());
3114
3115 result->Print();
3116
3117 Double_t value = result->GetValue(resultName.Data());
3118 Double_t error = result->GetErrorStat(resultName.Data());
3119
3120 g->SetPoint(i,runNumber,value);
3121 g->SetPointError(i,1,error);
3122
3123 Double_t n = mreal.CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL/run:%d",runNumber));
3124
3125 weightedMean += n*result->GetValue(resultName.Data());
3126 sumOfWeights += n;
3127
3128 mean += result->GetValue(resultName.Data());
3129
3130// std::cout << result->SubResults() << " " << result->GetError(resultName.Data()) << std::endl;
3131
2f331ac9 3132 }
3133
a58729a5 3134 gStyle->SetOptFit(1);
3135 g->Draw("alp");
3136 g->Fit("pol0");
3137 g->SetTitle("");
3138 g->GetXaxis()->SetTitle("Run number");
3139 g->GetXaxis()->SetNoExponent();
3140 if ( TString(what) == "Y" )
3141 {
3142 g->GetYaxis()->SetTitle("J/#psi yield");
3143 }
3144 else if ( TString(what) == "R" )
3145 {
3146 g->GetYaxis()->SetTitle(Form("R_{%s}^{J/#psi}",direction.Data()));
3147 }
3148 else if ( TString(what).Contains("Acc") )
3149 {
3150 g->GetYaxis()->SetTitle("Acc#timesEff_{J/#psi}");
3151 }
3152 else if ( TString(what).Contains("MBR") )
3153 {
3154 g->GetYaxis()->SetTitle("CINT7 / CINT7&0MUL");
3155 }
2f331ac9 3156
a58729a5 3157 if (CompactGraphs())
3158 {
3159 Compact(*g);
3160 }
2f331ac9 3161
a58729a5 3162 mean /= runs.size();
3163 weightedMean /= sumOfWeights;
2f331ac9 3164
a58729a5 3165 AliInfoClass(Form("Mean %e Weighted Mean %e",mean,weightedMean));
3166
3167 return g;
3168}
3169
3170//______________________________________________________________________________
3171void AliAnalysisMuMu::GetMBR(Int_t runNumber, const char* eventSelection, Double_t& value, Double_t& error) const
3172{
3173 // Get the scaling factor to go from CMUL to CINT7 for a given event selection
3174 value = 0.0;
3175 error = 0.0;
3176 if ( strlen(eventSelection) > 0 )
2f331ac9 3177 {
81dfe194 3178 ULong64_t a = TMath::Nint(fCounterCollection->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:%s/run:%d",
3179 eventSelection,runNumber)));
2f331ac9 3180
81dfe194 3181 ULong64_t b = TMath::Nint(fCounterCollection->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:%s/run:%d",
3182 eventSelection,runNumber)));
2f331ac9 3183
a58729a5 3184 value = b > 0 ? a/b : 0;
3185 error = value*AliAnalysisMuMuResult::ErrorAB(a,TMath::Sqrt(a),b,TMath::Sqrt(b));
2f331ac9 3186 }
a58729a5 3187}
3188
3189//______________________________________________________________________________
3190TGraph* AliAnalysisMuMu::MBREvolution(const char* eventSelection1, const char* eventSelection2) const
3191{
3192 if (!fCounterCollection) return 0x0;
2f331ac9 3193
a58729a5 3194 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
3195 runs->Sort();
3196 TIter nextRun(runs);
3197 TObjString* srun;
2f331ac9 3198
a58729a5 3199 std::set<int> runnumbers;
3200
3201 TGraphErrors* g = new TGraphErrors(runs->GetEntries());
3202 Int_t i(0);
2f331ac9 3203
a58729a5 3204 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2f331ac9 3205 {
a58729a5 3206 Int_t runNumber = srun->String().Atoi();
3207
3208 runnumbers.insert(runNumber);
3209
3210 Double_t mbr1,mbrError1;
3211 Double_t mbr2,mbrError2;
3212
3213 GetMBR(runNumber,eventSelection1,mbr1,mbrError1);
3214
3215 GetMBR(runNumber,eventSelection2,mbr2,mbrError2);
3216
3217 Double_t mbr = mbr1;
3218
3219 if ( mbr2 > 0 ) mbr /= mbr2;
3220
3221 Double_t mbrError = mbr*AliAnalysisMuMuResult::ErrorAB(mbr1,mbrError1,mbr2,mbrError2);
3222
3223 g->SetPoint(i,runNumber,mbr);
3224 g->SetPointError(i,0.5,mbrError);
3225
3226 ++i;
2f331ac9 3227 }
3228
a58729a5 3229 g->GetXaxis()->SetNoExponent();
3230
3231 AliAnalysisTriggerScalers ts(RunNumbers(),fgOCDBPath.Data());
3232 ts.DrawFills(0,1000);
2f331ac9 3233
a58729a5 3234 return g;
2f331ac9 3235}
a58729a5 3236
3237
3238
3239