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