1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 #include "AliAnalysisMuMu.h"
21 #include "AliAnalysisMuMuBinning.h"
22 #include "AliAnalysisMuMuFnorm.h"
23 #include "AliAnalysisMuMuGraphUtil.h"
24 #include "AliAnalysisMuMuJpsiResult.h"
25 #include "AliAnalysisMuMuSpectra.h"
26 #include "AliAnalysisTriggerScalers.h"
27 #include "AliCounterCollection.h"
28 #include "AliHistogramCollection.h"
30 #include "AliMergeableCollection.h"
31 #include "Riostream.h"
32 #include "TArrayL64.h"
39 #include "TGraphErrors.h"
43 #include "THashList.h"
46 #include "TLegendEntry.h"
51 #include "TObjArray.h"
52 #include "TObjString.h"
53 #include "TParameter.h"
54 #include "TPaveText.h"
56 #include "TStopwatch.h"
64 ClassImp(AliAnalysisMuMu)
66 TString AliAnalysisMuMu::fgOCDBPath("raw://");
68 TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-B-NOPF-MUON");
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");
72 TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON");
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");
76 TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD");
78 //,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
80 TString AliAnalysisMuMu::fgDefaultEventSelectionList("PSALL"); // for real data, for simulation see AliAnalysisMuMu ctor
82 TString AliAnalysisMuMu::fgDefaultPairSelectionList("pMATCHLOWRABSBOTH");
84 TString AliAnalysisMuMu::fgDefaultCentralitySelectionList("PP");
86 //TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOW:2,PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
87 TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
89 TString AliAnalysisMuMu::fgDefaultEventSelectionForSimulations("ALL");
90 TString AliAnalysisMuMu::fgDefaultDimuonTriggerForSimulations("CMULLO-B-NOPF-MUON");
92 Bool_t AliAnalysisMuMu::fgIsCompactGraphs(kFALSE);
94 //_____________________________________________________________________________
95 TString First(const TString s)
99 TObjArray* tokens = s.Tokenize(",");
101 if (!tokens) return rv;
103 rv = static_cast<TObjString*>(tokens->First())->String();
110 //_____________________________________________________________________________
111 TString FindTrigger(const AliMergeableCollection& mc,
113 const char* selection,
115 const char* centrality)
117 /// find the trigger containing the MinvPt histograms
119 std::vector<std::string> trigger2test;
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("-") )
126 trigger2test.push_back(base);
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));
135 trigger2test.push_back("ANY");
137 for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
139 std::string trigger = trigger2test[i];
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")
145 return trigger.c_str();
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 )
153 // AliWarningGeneral("FindTrigger",Form("tested trigger = %s",trigger2test[i].c_str()));
158 //_____________________________________________________________________________
159 AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName) : TObject(),
161 fCounterCollection(0x0),
162 fDimuonTriggers(fgDefaultDimuonTriggers),
163 fMuonTriggers(fgDefaultMuonTriggers),
164 fMinbiasTriggers(fgDefaultMinbiasTriggers),
165 fEventSelectionList(fgDefaultEventSelectionList),
166 fPairSelectionList(fgDefaultPairSelectionList),
167 fCentralitySelectionList(fgDefaultCentralitySelectionList),
168 fFitTypeList(fgDefaultFitTypeList),
170 fMergeableCollection(0x0),
172 fCorrectionPerRun(0x0),
173 fAssociatedSimulation(0x0)
177 GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
179 if ( fCounterCollection )
183 SetEventSelectionList("ALL");
184 SetDimuonTriggerList("CMULLO-B-NOPF-MUON");
185 SetFitTypeList("COUNTJPSI:1");
186 // SetFitTypeList("PSI1:1,COUNTJPSI:1");
189 if ( strlen(associatedSimFileName) )
191 fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
196 //_____________________________________________________________________________
197 AliAnalysisMuMu::~AliAnalysisMuMu()
201 if ( fAssociatedSimulation )
203 fAssociatedSimulation->Update();
208 delete fCounterCollection;
210 delete fMergeableCollection;
211 delete fCorrectionPerRun;
212 delete fAssociatedSimulation;
215 //_____________________________________________________________________________
216 void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
218 ULong64_t* totalNmsl,
219 ULong64_t* totalNmul)
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)
226 // filename is assumed to be a root filecontaining a list containing
227 // an AliCounterCollection (or directly an AliCounterCollection)
229 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
232 if (!fMergeableCollection || !fCounterCollection) return;
234 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
237 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
238 TIter nextTrigger(triggers);
240 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
242 Bool_t doPS = (events->FindObject("PSALL") != 0x0);
245 TObjString* strigger;
247 ULong64_t localNmb(0);
248 ULong64_t localNmsl(0);
249 ULong64_t localNmul(0);
251 if ( totalNmb) *totalNmb = 0;
252 if ( totalNmsl) *totalNmsl = 0;
253 if ( totalNmul ) *totalNmul = 0;
255 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
257 std::cout << Form("RUN %09d ",srun->String().Atoi());
268 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
271 if ( !fgDefaultMinbiasTriggers.Contains(strigger->String().Data()) &&
272 !fgDefaultMuonTriggers.Contains(strigger->String().Data()) &&
273 !fgDefaultDimuonTriggers.Contains(strigger->String().Data()) ) continue;
275 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
276 strigger->String().Data(),"ALL",srun->String().Atoi())));
278 details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
281 ULong64_t nps = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
282 strigger->String().Data(),"PSALL",srun->String().Atoi())));
286 details += TString::Format(" PS %5.1f %%",nps*100.0/n);
294 if ( fMinbiasTriggers.Contains(strigger->String()) )
297 if ( totalNmb) (*totalNmb) += n;
300 else if ( fMuonTriggers.Contains(strigger->String()) )
303 if ( totalNmsl) (*totalNmsl) += n;
306 else if ( fDimuonTriggers.Contains(strigger->String()) )
309 if ( totalNmul ) (*totalNmul) += n;
314 std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
315 nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
317 if ( detailTriggers )
319 std::cout << details.Data();
321 std::cout << std::endl;
324 if ( !totalNmul && !totalNmsl && !totalNmb )
326 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
327 localNmb,localNmsl,localNmul) << std::endl;
335 //_____________________________________________________________________________
336 void AliAnalysisMuMu::BasicCountsEvolution(const char* filelist, Bool_t detailTriggers)
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)
343 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
345 // To change the list of (single muon, dimuon, MB) triggers, use
346 // the SetDefault*TriggerList methods prior to call this one
349 TObjArray* files = ReadFileList(filelist);
351 if (!files || files->IsEmpty() ) return;
356 ULong64_t totalNmb(0);
357 ULong64_t totalNmsl(0);
358 ULong64_t totalNmul(0);
360 while ( ( str = static_cast<TObjString*>(next()) ) )
362 AliAnalysisMuMu m(str->String().Data());
368 m.BasicCounts(detailTriggers,&nmb,&nmsl,&nmul);
375 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
376 totalNmb,totalNmsl,totalNmul) << std::endl;
380 //_____________________________________________________________________________
381 void AliAnalysisMuMu::CentralityCheck(const char* filelist)
383 // Check if we get correctly filled centrality
385 TObjArray* files = ReadFileList(filelist);
387 if (!files || files->IsEmpty() ) return;
392 while ( ( str = static_cast<TObjString*>(next()) ) )
394 AliMergeableCollection* mc(0x0);
395 AliCounterCollection* cc(0x0);
396 AliAnalysisMuMuBinning* bin(0x0);
397 std::set<int> runnumbers;
399 if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
401 int run = RunNumberFromFileName(str->String().Data());
403 TH1* h = mc->Histo("/ALL/CPBI1MUL-B-NOPF-MUON/Centrality");
409 percent = 100*h->Integral(1,1)/h->Integral();
412 std::cout << Form("RUN %09d PERCENT %7.2f",run,percent) << std::endl;
421 //_____________________________________________________________________________
422 void AliAnalysisMuMu::CleanAllSpectra()
424 /// Delete all the spectra we may have
426 MC()->RemoveByType("AliAnalysisMuMuSpectra");
431 //_____________________________________________________________________________
432 TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
433 const char* backgroundresults)
435 TFile* fjpsi = FileOpen(jpsiresults);
436 TFile* fbck = FileOpen(backgroundresults);
438 if (!fjpsi || !fbck) return 0x0;
440 TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
442 std::vector<std::string> checks;
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");
449 if (!gjpsi) return 0x0;
451 TObjArray* a = new TObjArray;
454 for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
457 TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
459 if (!gback) continue;
461 if ( gjpsi->GetN() != gback->GetN() )
463 AliErrorClass("graphs have different number of points !");
467 TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
469 for ( int i = 0; i < gjpsi->GetN(); ++i )
473 gjpsi->GetPoint(i,r1,y1);
474 gback->GetPoint(i,r2,y2);
478 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
482 g->SetPoint(i,y2,y1);
483 // g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
486 g->SetMarkerStyle(25+j);
487 g->SetMarkerSize(1.2);
496 g->SetLineColor(j+1);
497 g->SetMarkerColor(j+1);
498 g->SetName(checks[j].c_str());
505 //_____________________________________________________________________________
506 TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
507 const char* simjpsiresults)
509 TFile* freal = FileOpen(realjpsiresults);
510 TFile* fsim = FileOpen(simjpsiresults);
512 if (!freal || !fsim) return 0x0;
514 TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
515 TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
517 TObjArray* a = new TObjArray;
520 if ( greal->GetN() != gsim->GetN() )
522 AliErrorClass("graphs have different number of points !");
526 TGraphErrors* g = new TGraphErrors(greal->GetN());
527 TGraphErrors* gratio = new TGraphErrors(greal->GetN());
529 for ( int i = 0; i < greal->GetN(); ++i )
533 greal->GetPoint(i,r1,y1);
534 gsim->GetPoint(i,r2,y2);
538 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
544 if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
547 g->SetPointError(i,0,0);
551 g->SetPoint(i,y2,y1);
552 g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
555 gratio->SetPoint(i,r1,ratio);
558 g->SetMarkerStyle(25);
559 g->SetMarkerSize(1.2);
566 g->SetMarkerColor(1);
567 g->SetName("jpsipercmuurealvssim");
572 gsim->SetLineColor(4);
582 //_____________________________________________________________________________
583 TObjArray* AliAnalysisMuMu::ComputeBackgroundEvolution(const char* filelist,
584 const char* triggerList,
586 const char* outputFile,
587 const char* outputMode)
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.
593 TObjArray* files = ReadFileList(filelist);
595 if (!files || files->IsEmpty() ) return 0x0;
600 const char* ps = "PSALL";
601 const char* centrality = "PP";
602 const char* ts1 = "sMATCHLOWRABS";
603 const char* ts2 = "sMATCHLOWRABSDCA";
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;
610 TObjArray* triggers = TString(triggerList).Tokenize(",");
612 TObjArray* a = new TObjArray;
615 Bool_t bothSigns(kFALSE);
617 while ( ( str = static_cast<TObjString*>(next()) ) )
619 AliInfoClass(str->String().Data());
621 AliMergeableCollection* mc(0x0);
622 AliCounterCollection* cc(0x0);
623 AliAnalysisMuMuBinning* bin(0x0);
624 std::set<int> runnumbers;
626 if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
628 TIter nextObject(mc->CreateIterator());
630 int nplus(0), nminus(0);
632 while ( ( o = nextObject() ) )
634 if ( o->InheritsFrom("TH1") )
639 TH1* h = static_cast<TH1*>(o);
641 if ( TString(h->GetName()).EndsWith("Plus") )
645 if ( TString(h->GetName()).EndsWith("Minus") )
651 if (nminus==nplus && nplus>0 )
656 AliInfoClass(Form("Both signs = %d",bothSigns));
658 TIter nextTrigger(triggers);
666 while ( ( trigger = static_cast<TObjString*>(nextTrigger()) ) )
670 format = "/%s/%s/%s/%s/PtEtaMuPlus:py";
674 format = "/%s/%s/%s/%s/PtEtaMu:py";
677 TString hname(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
679 h1p = mc->Histo(hname.Data());
683 AliInfoClass(Form("Could not get %s",hname.Data()));
687 AliInfoClass(Form("Will use trigger %s",trigger->String().Data()));
689 h2p = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
693 format.ReplaceAll("Plus","Minus");
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));
703 if (h1m && h2m && h1p && h2p)
705 Int_t bin1 = h1m->GetXaxis()->FindBin(ptmin);
706 Int_t bin2 = h1m->GetXaxis()->GetNbins();
708 runs[trigger->String().Data()].push_back(RunNumberFromFileName(str->String().Data()));
709 errruns[trigger->String().Data()].push_back(0.5);
712 double v1 = h2m->IntegralAndError(bin1,bin2,e1);
713 double v2 = h1m->IntegralAndError(bin1,bin2,e2);
714 double value = 100*(1.0-v1/v2);
717 yminus[trigger->String().Data()].push_back(value);
718 // double e1 = 1.0/TMath::Sqrt(h1m->GetEntries());
719 // double e2 = 1.0/TMath::Sqrt(h2m->GetEntries());
720 erryminus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
722 v1=h2p->IntegralAndError(bin1,bin2,e1);
723 v2=h1p->IntegralAndError(bin1,bin2,e1);
724 value = 100*(1.0-v1/v2);
727 yplus[trigger->String().Data()].push_back(value);
728 // e1 = 1.0/TMath::Sqrt(h1p->GetEntries());
729 // e2 = 1.0/TMath::Sqrt(h2p->GetEntries());
730 erryplus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
734 std::cout << Form("Error : h1m %p h2m %p h1p %p h2p %p",h1m,h2m,h1p,h2p) << std::endl;
740 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(str->String().Data()));
746 TFile* f = new TFile(outputFile,outputMode);
748 std::map<std::string, std::vector<float> >::const_iterator it;
750 for ( it = runs.begin(); it != runs.end(); ++it )
752 std::string triggerName = it->first;
754 TGraphErrors* gp = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yplus[triggerName][0],&errruns[triggerName][0],&erryplus[triggerName][0]);
755 TGraphErrors* gm(0x0);
759 gm = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yminus[triggerName][0],&errruns[triggerName][0],&erryminus[triggerName][0]);
764 gp->Write(Form("muplus_%s",triggerName.c_str()),TObject::kOverwrite);
765 gm->Write(Form("muminus_%s",triggerName.c_str()),TObject::kOverwrite);
769 gp->Write(Form("mu_%s",triggerName.c_str()),TObject::kOverwrite);
779 //_____________________________________________________________________________
781 AliAnalysisMuMu::ComputeJpsiEvolution(const char* filelist, const char* triggerList,
782 const char* outputFile)
784 /// Compute some jpsi information for a list of files / trigger combinations
786 TObjArray* files = ReadFileList(filelist);
788 if (!files || files->IsEmpty() ) return 0x0;
790 TMap results; // one TObjString->TObjArray per file
791 results.SetOwnerKeyValue(kTRUE,kTRUE);
793 TIter nextFile(files);
797 // while ( ( str = static_cast<TObjString*>(nextFile()) ) )
799 // std::cout << str->String().Data() << std::endl;
801 // AliAnalysisMuMu m(str->String().Data());
803 // m.SetDimuonTriggerList(triggerList);
805 // TMap* map = m.Jpsi();
809 // AliWarningClass(Form("Got no jpsi for %s",str->String().Data()));
812 // results.Add(new TObjString(str->String()), map);
815 // if (!results.GetSize()) return 0x0;
817 // compute the total over all files
819 TMap* total = new TMap;
820 total->SetOwnerKeyValue(kTRUE,kTRUE);
823 TObjArray* triggers = TString(triggerList).Tokenize(",");
825 TIter nextTrigger(triggers);
826 TObjString* trigger(0x0);
828 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
833 AliAnalysisMuMuResult* ref(0x0);
835 while ( ( str = static_cast<TObjString*>(nextFile()) ) )
837 // TObjArray* a = static_cast<TObjArray*>(results.GetValue(str->String().Data()));
839 AliInfoClass("FIXME: write the merging of AliAnalysisMuMuResult objects !");
841 // AliAnalysisMuMuResult* r(0x0);
845 // r = static_cast<AliAnalysisMuMuResult*>(a->FindObject(trigger->String().Data()));
849 // if (!ref) ref = r;
853 // TH1* htmp = static_cast<TH1*>(r->Minv()->Clone());
865 // n += r->NofTriggers();
876 // AliAnalysisMuMuResult* sum = new AliAnalysisMuMuResult(*hminv,
877 // ref->TriggerClass(),
878 // ref->EventSelection(),
879 // ref->PairSelection(),
880 // ref->CentralitySelection(),
881 // AliAnalysisMuMuBinning::Range());
883 // sum->SetNofTriggers(n);
885 // sum->SetNofRuns(nruns);
889 // total->Add(new TObjString(trigger->String().Data()),sum);
893 AliInfoClass("--------------------------------------");
894 StdoutToAliInfoClass(total->Print(););
896 AliInfoClass("---------------------------Going to write file");
898 TFile* fout = new TFile(outputFile,"RECREATE");
900 results.Write("rbr",TObject::kSingleKey);
902 total->Write("total",TObject::kSingleKey);
908 AliInfoClass(Form("%d files analyzed",files->GetEntries()));
913 //_____________________________________________________________________________
914 Bool_t AliAnalysisMuMu::DecodeFileName(const char* filename,
920 // tries to extract period and pass numbers from a file name.
922 esdpass=aodtrain=runnumber=-1;
925 TString sfile(gSystem->BaseName(filename));
927 if (!sfile.BeginsWith("LHC") && !sfile.BeginsWith("SIM") )
929 // does not look nice but let's try at least to get the period
931 TObjArray* tmp = sfile.Tokenize(".");
933 for ( Int_t j = 0; j < tmp->GetEntries(); ++j )
935 TString s = static_cast<TObjString*>(tmp->At(j))->String();
937 if ( s.BeginsWith("lhc") )
940 period.ReplaceAll("lhc","LHC");
947 if ( period.Length() )
952 std::cerr << Form("filename %s does not start with LHC or SIM",filename) << std::endl;
960 if ( sfile.BeginsWith("LHC") )
962 if ( sfile.Contains("pass") && sfile.Contains("AODMUON") )
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);
969 else if ( sfile.Contains("pass") && sfile.Contains("_muon_") && sfile.Contains("AOD000") )
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);
976 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AODMUON000") )
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);
983 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AOD000") )
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);
990 else if ( sfile.Contains("AODMUON" ) )
992 sscanf(sfile.Data(),"LHC%2d%c_AODMUON%03d_%09d",&year,&p,&aodtrain,&runnumber);
993 period = TString::Format("LHC%2d%c",year,p);
996 else if ( sfile.Contains("AOD000") )
998 sscanf(sfile.Data(),"LHC%2d%c_muons_AOD000_%09d",&year,&p,&runnumber);
999 period = TString::Format("LHC%2d%c",year,p);
1002 else if ( sfile.Contains("ESD_OUTER000"))
1004 sscanf(sfile.Data(),"LHC%2d%c_cpass1_ESD_OUTER000_%09d",&year,&p,&runnumber);
1008 else if ( sfile.BeginsWith("SIM_JPSI3" ) )
1010 sscanf(sfile.Data(),"SIM_JPSI3_%09d",&runnumber);
1013 else if ( sfile.BeginsWith("SIM_UPSILON" ) )
1015 sscanf(sfile.Data(),"SIM_UPSILON_%09d",&runnumber);
1018 else if ( sfile.BeginsWith("SIM_JPSI" ) )
1020 sscanf(sfile.Data(),"SIM_JPSI_LHC%2d%c_%09d",&year,&p,&runnumber);
1021 period = TString::Format("LHC%2d%c",year,p);
1027 std::cerr << Form("Can not decode %s",filename) << std::endl;
1034 //_____________________________________________________________________________
1035 void AliAnalysisMuMu::DrawMinv(const char* type,
1036 const char* particle,
1037 const char* trigger,
1038 const char* eventType,
1039 const char* pairCut,
1040 const char* centrality,
1041 const char* subresultname,
1042 const char* flavour) const
1044 /// Draw minv spectra for binning of given type
1046 if (!MC() || !BIN()) return;
1048 TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
1051 AliError(Form("Could not get %s bins",type));
1058 TString sparticle(particle);
1059 if ( sparticle=="PSI" )
1068 Int_t n = bins->GetEntries();
1076 ny = TMath::Nint(TMath::Sqrt(n));
1080 TString stype(type);
1083 TString spectraName(Form("/%s/%s/%s/%s/%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data()));
1085 if ( strlen(flavour))
1088 spectraName += flavour;
1091 AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(spectraName.Data()));
1093 AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
1095 TObjArray* spectraBins(0x0);
1098 spectraBins = spectra->BinContentArray();
1101 TCanvas* c = new TCanvas;
1104 gStyle->SetOptFit(1112);
1106 c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
1107 (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
1111 AliAnalysisMuMuBinning::Range* r;
1114 while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1116 TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
1118 AliDebug(1,name.Data());
1120 AliAnalysisMuMuJpsiResult* spectraBin(0x0);
1124 AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
1126 spectraBin = static_cast<AliAnalysisMuMuJpsiResult*>(sr->SubResult(subresultname));
1128 AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
1131 TH1* h = MC()->Histo(name.Data());
1135 h = spectraBin->Minv();
1145 h->GetXaxis()->SetRangeUser(xmin,xmax);
1149 TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
1158 TObject* stats = h->FindObject("stats");
1161 stats->Draw("same");
1169 //_____________________________________________________________________________
1170 void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
1172 /// Draw minv spectra for binning of given type
1174 DrawMinv(type,particle,
1175 First(DimuonTriggerList()).Data(),
1176 First(EventSelectionList()).Data(),
1177 First(PairSelectionList()).Data(),
1178 First(CentralitySelectionList()).Data(),
1183 //___________________________________________________________________
1184 void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
1186 // Actions in reponse to mouse button events.
1188 TCanvas* c = static_cast<TCanvas*>(gTQSender);
1189 TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
1192 // if ((event == kButton1Down) ||
1193 if (event == kButton1Double)
1196 // Float_t x = pad->AbsPixeltoX(px);
1197 // Float_t y = pad->AbsPixeltoY(py);
1198 // x = pad->PadtoX(x);
1199 // y = pad->PadtoY(y);
1201 // std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
1203 if ( sel && sel->InheritsFrom("TH1") )
1205 TCanvas* clocal = new TCanvas;
1212 TList* list = pad->GetListOfPrimitives();
1216 while ( ( h = next() ) )
1218 if ( h->InheritsFrom("TH1") )
1220 TCanvas* clocal = new TCanvas;
1230 // std::cout << std::endl;
1237 //_____________________________________________________________________________
1239 AliAnalysisMuMu::ExpandPathName(const char* file)
1241 // An expand method that lives alien URL as they are
1244 if ( !sfile.BeginsWith("alien://") )
1246 return gSystem->ExpandPathName(file);
1250 if (!gGrid) TGrid::Connect("alien://");
1251 if (!gGrid) return "";
1257 //_____________________________________________________________________________
1258 void AliAnalysisMuMu::TwikiOutputFnorm(const char* series) const
1260 // make a twiki-compatible output of the Fnorm factor(s)
1261 TObjArray* what = TString(series).Tokenize(",");
1266 std::cout << "| *Run* |";
1267 while ( ( s = static_cast<TObjString*>(next())) )
1269 TGraph* g = static_cast<TGraph*>(MC()->GetObject(Form("/FNORM/GRAPHS/%s",s->String().Data())));
1272 AliError(Form("Could not find graph for %s",s->String().Data()));
1275 std::cout << " *" << s->String().Data();
1276 if ( s->String().BeginsWith("RelDif") ) std::cout << " %";
1281 std::cout << std::endl;
1283 TGraphErrors* g0 = static_cast<TGraphErrors*>(graphs.First());
1286 for ( Int_t i = 0; i < g0->GetN(); ++i )
1290 msg.Form("|%6d|",TMath::Nint(g0->GetX()[i]));
1292 for ( Int_t j = 0; j < graphs.GetEntries(); ++j )
1294 TGraphErrors* g = static_cast<TGraphErrors*>(graphs.At(j));
1296 msg += TString::Format(" %6.2f +- %6.2f |",g->GetY()[i],g->GetEY()[i]);
1299 std::cout << msg.Data() << std::endl;
1304 std::cout << "|*Weigthed mean (*)*|";
1306 AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1310 AliError("Could not find Fnorm result !");
1315 while ( ( s = static_cast<TObjString*>(next())) )
1317 TString var("Fnorm");
1320 if ( s->String().BeginsWith("Fnorm") )
1322 r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1324 else if ( s->String().BeginsWith("RelDif") )
1326 r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/RelDif"));
1331 r->Include(s->String().Data());
1333 std::cout << Form(" * %5.2f +- %5.2f %s * |",
1334 r->GetValue(var.Data()),
1335 r->GetErrorStat(var.Data()),
1341 std::cout << std::endl;
1343 std::cout << "|*RMS*|";
1345 while ( ( s = static_cast<TObjString*>(next())) )
1347 TString var("Fnorm");
1349 if ( s->String().BeginsWith("Fnorm") )
1351 r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1353 else if ( s->String().BeginsWith("RelDif") )
1355 r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/RelDif"));
1359 r->Include(s->String().Data());
1361 Double_t d = 100.0*r->GetRMS(var.Data())/r->GetValue(var.Data());
1363 std::cout << Form(" * %5.2f (%5.2f %%) * |",
1364 r->GetRMS(var.Data()),d);
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;
1373 //_____________________________________________________________________________
1374 void AliAnalysisMuMu::FigureOutputFnorm(const char* filelist)
1376 /// Make some figure of the Fnorm factors for files in filelist
1378 TObjArray* periods = ReadFileList(filelist);
1380 if (!periods || periods->IsEmpty() ) return;
1382 TIter next(periods);
1384 TObjArray fnormoffline1;
1385 TObjArray fnormoffline2;
1387 TObjArray correctionPSMUL;
1388 TObjArray correctionPSMB;
1389 TObjArray correctionPUPS;
1390 TObjArray correctionPSRatio;
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);
1400 for ( Int_t i = 0; i <= periods->GetLast(); ++i )
1402 TString period("unknown");
1404 TString filename(static_cast<TObjString*>(periods->At(i))->String());
1408 if (!DecodeFileName(filename,period,dummy,dummy,dummy))
1413 if ( gSystem->AccessPathName(filename) )
1415 AliErrorClass(Form("Could not find file %s. Skipping it.",filename.Data()));
1420 AliAnalysisMuMu m(filename.Data());
1422 fnormoffline1.Add(m.MC()->GetObject("/FNORM/GRAPHS/FnormOffline1PUPS")->Clone());
1423 fnormoffline2.Add(m.MC()->GetObject("/FNORM/GRAPHS/FnormOffline2PUPS")->Clone());
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());
1430 correctionPSRatio.Add(m.MC()->GetObject("/FNORM/GRAPHS/CorrectionPSRatio")->Clone());
1433 reldif.Add(m.MC()->GetObject("/FNORM/GRAPHS/RelDifFnormScalersPUPSvsFnormOffline2PUPS")->Clone());
1437 AliAnalysisMuMuGraphUtil gu(fgOCDBPath);
1439 gu.ShouldDrawPeriods(kTRUE);
1443 a.Add(gu.Combine(fnormoffline1,fgIsCompactGraphs));
1444 a.Add(gu.Combine(fnormoffline2,fgIsCompactGraphs));
1447 Double_t ymax(3000.0);
1449 new TCanvas("fnormoffline","fnormoffline");
1451 gu.PlotSameWithLegend(a,ymin,ymax);
1453 new TCanvas("corrections","corrections");
1457 a.Add(gu.Combine(correctionPSMB,fgIsCompactGraphs));
1458 a.Add(gu.Combine(correctionPSMUL,fgIsCompactGraphs));
1459 a.Add(gu.Combine(correctionPUPS,fgIsCompactGraphs));
1461 gu.PlotSameWithLegend(a,0.8,1.2);
1463 new TCanvas("psratio","psratio");
1467 a.Add(gu.Combine(correctionPSRatio,fgIsCompactGraphs));
1469 gu.PlotSameWithLegend(a,0.8,1.2);
1471 new TCanvas("offlinevsscalers","fig:offlinevsscalers");
1475 a.Add(gu.Combine(reldif,fgIsCompactGraphs));
1477 gu.PlotSameWithLegend(a,-15,15);
1480 //_____________________________________________________________________________
1482 AliAnalysisMuMu::FileOpen(const char* file)
1484 // Open a file after expansion of its name
1486 return TFile::Open(ExpandPathName(file).Data());
1489 //_____________________________________________________________________________
1490 TString AliAnalysisMuMu::First(const TString& list) const
1492 TObjArray* a = list.Tokenize(",");
1493 if ( a->GetLast() < 0 ) return "";
1495 TString rv = static_cast<TObjString*>(a->First())->String();
1502 //_____________________________________________________________________________
1503 AliAnalysisMuMuSpectra*
1504 AliAnalysisMuMu::FitParticle(const char* particle,
1505 const char* trigger,
1506 const char* eventType,
1507 const char* pairCut,
1508 const char* centrality,
1509 const AliAnalysisMuMuBinning& binning)
1511 // Fit the minv spectra to find the given particle
1512 // Returns an array of AliAnalysisMuMuResult objects
1516 TObjArray* bins = binning.CreateBinObjArray(particle);
1519 AliError(Form("Did not get any bin for particle %s",particle));
1523 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1524 if ( !triggers->FindObject(trigger) )
1526 AliDebug(1,Form("Did not find trigger %s",trigger));
1533 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1534 if ( !events->FindObject(eventType) )
1536 AliError(Form("Did not find eventType %s",eventType));
1543 Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
1547 AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
1552 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
1553 Int_t nruns = runs->GetEntries();
1559 AliAnalysisMuMuSpectra* spectra(0x0);
1561 AliAnalysisMuMuBinning::Range* bin;
1564 TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
1565 TIter nextFitType(fitTypeArray);
1566 TObjString* fitType;
1569 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1571 TString hname(Form("MinvUS%s",bin->AsString().Data()));
1573 TH1* hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),hname.Data());
1577 if (!fBinning && bin->IsNullObject() )
1579 // old file, we only had MinvUSPt
1580 hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
1585 AliDebug(1,Form("Could not find histo %s",hname.Data()));
1590 hminv = static_cast<TH1*>(hminv->Clone(Form("minv%d",n++)));
1593 AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(*hminv,
1600 r->SetNofTriggers(ntrigger);
1601 r->SetNofRuns(nruns);
1603 nextFitType.Reset();
1605 while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
1607 AliDebug(1,Form("<<<<<< fitType=%s bin=%s",fitType->String().Data(),bin->Flavour().Data()));
1609 if ( fitType->String().BeginsWith("PSILOWMCTAILS") )
1611 std::vector<Double_t> par;
1612 par = GetMCCB2Tails(*bin);
1615 r->AddFit(fitType->String().Data(),par.size(),&par[0]);
1620 r->AddFit(fitType->String());
1624 flavour = bin->Flavour();
1628 TString spectraName(binning.GetName());
1629 if ( flavour.Length() > 0 )
1632 spectraName += flavour;
1634 spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
1637 spectra->AdoptResult(*bin,r);
1639 if ( IsSimulation() )
1641 SetNofInputParticles(*r);
1647 delete fitTypeArray;
1653 //_____________________________________________________________________________
1654 std::vector<Double_t>
1655 AliAnalysisMuMu::GetMCCB2Tails(const AliAnalysisMuMuBinning::Range& bin) const
1657 /// Get the tails from the associated simulation
1659 std::vector<Double_t> par;
1663 AliError("Cannot get MC tails without an associated simulation file !");
1668 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(SIM()->GetSpectra(bin.Type().Data(),bin.Flavour().Data()));
1672 AliError(Form("Could not find spectra %s,%s for associated simulation",bin.Type().Data(),bin.Flavour().Data()));
1673 fAssociatedSimulation->MC()->Print("*:Ali*");
1678 AliDebug(1,Form("AliAnalysisMuMuSpectra* s = reinterpret_cast<AliAnalysisMuMuSpectra*>(%p)",s));
1681 AliAnalysisMuMuResult* r = s->GetResultForBin(bin);
1685 AliAnalysisMuMuJpsiResult* r1 = dynamic_cast<AliAnalysisMuMuJpsiResult*>(r->SubResult("JPSI:1"));
1688 TF1* func = static_cast<TF1*>(r1->Minv()->GetListOfFunctions()->FindObject("fitTotal"));
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"));
1702 //_____________________________________________________________________________
1703 ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
1705 // Get the expected (from OCDB scalers) trigger count
1707 AliAnalysisTriggerScalers ts(runNumber,fgOCDBPath.Data());
1709 TObjArray* triggers = TString(triggerList).Tokenize(",");
1710 TObjString* trigger;
1711 TIter next(triggers);
1714 while ( ( trigger = static_cast<TObjString*>(next()) ) )
1716 AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
1728 //_____________________________________________________________________________
1729 AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
1731 /// get a given spectra
1733 TString swhat(what);
1734 TString sflavour(flavour);
1738 TString spectraName(Form("/PSALL/%s/PP/%s/PSI-%s",
1739 First(fDimuonTriggers).Data(),
1740 First(fPairSelectionList).Data(),
1743 if (sflavour.Length()>0)
1746 spectraName += sflavour.Data();
1751 spectraName.ReplaceAll("PSALL",fgDefaultEventSelectionForSimulations.Data());
1752 spectraName.ReplaceAll(First(fgDefaultDimuonTriggers).Data(),fgDefaultDimuonTriggerForSimulations.Data());
1755 return SPECTRA(spectraName.Data());
1758 //_____________________________________________________________________________
1759 UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
1760 const char* eventSelection, Int_t runNumber)
1762 TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
1763 TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
1764 TObjArray* a = TString(triggerList).Tokenize(" ");
1770 TString sEventSelection(eventSelection);
1771 sEventSelection.ToUpper();
1773 if ( kevent->FindObject(sEventSelection.Data()) )
1775 while ( ( str = static_cast<TObjString*>(next()) ) )
1777 if ( ktrigger->FindObject(str->String().Data()) )
1779 if ( runNumber < 0 )
1781 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
1785 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
1797 //_____________________________________________________________________________
1799 AliAnalysisMuMu::GetCollections(const char* rootfile,
1800 AliMergeableCollection*& mc,
1801 AliCounterCollection*& cc,
1802 AliAnalysisMuMuBinning*& bin,
1803 std::set<int>& runnumbers)
1809 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
1812 f = TFile::Open(rootfile);
1815 if ( !f || f->IsZombie() )
1820 f->GetObject("MC",mc);
1821 f->GetObject("CC",cc);
1823 TIter next(f->GetListOfKeys());
1826 while ( ( key = static_cast<TKey*>(next())) && !bin )
1828 if ( strcmp(key->GetClassName(),"AliAnalysisMuMuBinning")==0 )
1830 bin = dynamic_cast<AliAnalysisMuMuBinning*>(key->ReadObj());
1836 AliErrorClass("Old file. Please upgrade it!");
1842 TObjArray* runs = cc->GetKeyWords("run").Tokenize(",");
1844 TIter nextRun(runs);
1849 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
1851 runnumbers.insert(srun->String().Atoi());
1859 //_____________________________________________________________________________
1860 Bool_t AliAnalysisMuMu::IsSimulation() const
1862 // whether or not we have MC information
1864 if (!fMergeableCollection) return kFALSE;
1866 return ( fMergeableCollection->Histo(Form("/INPUT/%s/MinvUS",fgDefaultEventSelectionForSimulations.Data())) != 0x0 );
1869 //_____________________________________________________________________________
1871 AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour)
1873 // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
1874 // what="integrated" => fit only fully integrated MinvUS
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
1881 if (!fMergeableCollection)
1883 AliError("No mergeable collection. Consider Upgrade()");
1889 TObjArray* triggerArray = fDimuonTriggers.Tokenize(",");
1890 TObjArray* eventTypeArray = fEventSelectionList.Tokenize(",");
1891 TObjArray* pairCutArray = fPairSelectionList.Tokenize(",");
1892 TObjArray* whatArray = TString(what).Tokenize(",");
1894 TIter nextTrigger(triggerArray);
1895 TIter nextEventType(eventTypeArray);
1896 TIter nextPairCut(pairCutArray);
1897 TIter nextWhat(whatArray);
1899 TObjString* trigger;
1900 TObjString* eventType;
1901 TObjString* pairCut;
1904 while ( ( swhat = static_cast<TObjString*>(nextWhat()) ) )
1906 AliAnalysisMuMuBinning* binning(0x0);
1908 if ( fBinning && swhat->String().Length() > 0 )
1910 binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
1914 binning = new AliAnalysisMuMuBinning;
1915 binning->AddBin("psi",swhat->String().Data());
1918 StdoutToAliDebug(1,std::cout << "++++++++++++ swhat=" << swhat->String().Data() << std::endl;);
1922 AliError("oups. binning is NULL");
1926 StdoutToAliDebug(1,binning->Print(););
1928 nextTrigger.Reset();
1930 while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
1932 AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
1934 nextEventType.Reset();
1936 while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
1938 AliDebug(1,Form("--EVENTTYPE %s",eventType->String().Data()));
1940 nextPairCut.Reset();
1942 while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
1944 AliDebug(1,Form("----PAIRCUT %s",pairCut->String().Data()));
1946 AliDebug(1,"----Fitting...");
1948 AliAnalysisMuMuSpectra* spectra = FitParticle("psi",
1949 trigger->String().Data(),
1950 eventType->String().Data(),
1951 pairCut->String().Data(),
1955 AliDebug(1,Form("----fitting done spectra = %p",spectra));
1961 TString id(Form("/%s/%s/PP/%s",eventType->String().Data(),
1962 trigger->String().Data(),
1963 pairCut->String().Data()));
1965 TObject* o = fMergeableCollection->GetObject(id.Data(),spectra->GetName());
1967 AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
1971 AliWarning(Form("Replacing %s/%s",id.Data(),spectra->GetName()));
1972 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectra->GetName()));
1975 fMergeableCollection->Adopt(id.Data(),spectra);
1977 StdoutToAliDebug(1,spectra->Print(););
1985 delete triggerArray;
1986 delete eventTypeArray;
1987 delete pairCutArray;
1989 StdoutToAliDebug(1,timer.Print(););
1994 // ReOpen(fFilename,"UPDATE");
1995 // fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
1996 // ReOpen(fFilename,"READ");
2004 //_____________________________________________________________________________
2005 void AliAnalysisMuMu::LatexOutputFnorm(const char* filelist, const char* subresultnames, Bool_t rms)
2007 /// Make a LaTeX output of the Fnorm factors for each file in filelist
2009 TObjArray* periods = ReadFileList(filelist);
2011 if (!periods || periods->IsEmpty() ) return;
2013 TIter next(periods);
2015 TObjArray* subresults = TString(subresultnames).Tokenize(",");
2016 TIter nextSub(subresults);
2019 Int_t iclast = subresults->GetLast();
2022 std::cout << "\\begin{tabular}{l|r|r|r|r}" << std::endl;
2024 std::cout << "Period & ";
2029 while ( ( rname = static_cast<TObjString*>(nextSub())) )
2031 std::cout << rname->String().Data();
2040 std::cout << "\\\\" << std::endl << "\\hline" << std::endl;
2042 for ( Int_t i = 0; i <= periods->GetLast(); ++i )
2044 TString period("unknown");
2046 TString filename(static_cast<TObjString*>(periods->At(i))->String());
2050 if (!DecodeFileName(filename,period,dummy,dummy,dummy))
2055 if ( gSystem->AccessPathName(filename) )
2057 AliErrorClass(Form("Could not find file %s. Skipping it.",filename.Data()));
2062 AliAnalysisMuMu m(filename.Data());
2064 AliMergeableCollection* norm = m.MC()->Project("/FNORM/RESULTS/");
2066 AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Fnorm")->Clone());
2072 AliErrorClass(Form("Could not get result for file %s.",filename.Data()));
2076 std::cout << period << " & ";
2082 while ( ( rname = static_cast<TObjString*>(nextSub())) )
2084 TString name(rname->String());
2085 TString var("Fnorm");
2087 if ( name.BeginsWith("RelDif"))
2089 r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("RelDif"));
2092 else if ( name.BeginsWith("Correction"))
2094 r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Correction"));
2098 r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Fnorm"));
2101 r->Include(name.Data());
2105 std::cout << Form(" $%7.2f (%5.2f \%%)$ ",
2106 r->GetRMS(var.Data()),
2107 100.0*r->GetRMS(var.Data())/r->GetValue(var.Data()));
2112 std::cout << Form(" $%7.2f \\pm %5.3f$ ",r->GetValue(var.Data()),
2113 r->GetErrorStat(var.Data()));
2124 if ( i != periods->GetLast() )
2126 std::cout << "\\\\";
2129 std::cout << std::endl;
2135 std::cout << "\\end{tabular}" << std::endl;
2142 //_____________________________________________________________________________
2143 void AliAnalysisMuMu::PlotBackgroundEvolution(const char* gfile, const char* triggerList, Double_t ymax, Bool_t fillBoundaries)
2145 // plot the graphs found in the file (and generated using the ComputeBackgroundEvolution() method)
2147 TFile* f = TFile::Open(ExpandPathName(gfile).Data());
2149 if ( !f || !f->IsOpen() )
2157 TCanvas* c = new TCanvas("background-evolution","background-evolution");
2161 TLegend* l = new TLegend(0.4,0.6,0.97,0.97);
2163 l->SetTextColor(AliAnalysisMuMu::kBlue);
2164 l->SetLineColor(AliAnalysisMuMu::kBlue);
2166 TObjArray* triggers = TString(triggerList).Tokenize(",");
2168 gStyle->SetOptTitle(0);
2170 TObjString* str(0x0);
2171 TIter next(triggers);
2173 Int_t run1(99999999);
2176 std::set<int> runnumbers;
2178 while ( ( str = static_cast<TObjString*>(next()) ) )
2180 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2182 for ( Int_t ir = 0; ir < g->GetN(); ++ir )
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);
2191 AliInfoClass(Form("run1 %d run2 %d",run1,run2));
2195 TH2* hframe = new TH2F("hframe","hframe",run2-run1+1,run1,run2,100,ymin,ymax);
2197 hframe->GetXaxis()->SetNoExponent();
2198 hframe->GetYaxis()->SetTitle("Background percentage");
2202 AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath.Data());
2203 ts.DrawFills(ymin,ymax);
2208 while ( ( str = static_cast<TObjString*>(next()) ) )
2210 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2213 AliErrorClass(Form("Graph mu_%s not found",str->String().Data()));
2219 if (i==0) color = AliAnalysisMuMu::kBlue;
2220 if (i==1) color = AliAnalysisMuMu::kOrange;
2222 g->SetLineColor(color);
2223 g->SetMarkerColor(color);
2224 g->SetMarkerStyle(20+i);
2228 TLegendEntry* le = l->AddEntry(g,str->String().Data(),"lp");
2229 le->SetTextColor(color);
2231 g->GetYaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2232 g->GetXaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2238 hframe->Draw("sameaxis");
2243 //_____________________________________________________________________________
2245 AliAnalysisMuMu::PlotJpsiEvolution(const char* resultFile, const char* triggerList, Bool_t fillBoundaries,
2246 const char* efficiencyFile, Bool_t simulation)
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)
2251 std::map<int, float> efficiencies;
2253 if ( efficiencyFile && strlen(efficiencyFile) > 0 )
2255 std::ifstream in(gSystem->ExpandPathName(efficiencyFile));
2261 float dummy,errorl,errorh;
2263 while ( in.getline(line,1023,'\n') )
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);
2268 AliDebugClass(1,Form("%09d %8.6f +- %8.6f ",run,eff,errorl+errorh));
2270 efficiencies[run] = eff;
2275 TFile* f = TFile::Open(gSystem->ExpandPathName(resultFile));
2277 std::set<int> runnumbers;
2279 TMap* m = static_cast<TMap*>(f->Get("rbr"));
2285 files.SetOwner(kTRUE);
2287 while ( ( str = static_cast<TObjString*>(next())) )
2289 files.Add(new TObjString(str->String()));
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;
2299 TIter nextTrigger(TString(triggerList).Tokenize(","));
2300 TObjString* trigger(0x0);
2302 int runMin(100000000);
2305 TIter nextFile(&files);
2307 Double_t ymin(TMath::Limits<double>::Max());
2308 Double_t ymax(TMath::Limits<double>::Min());
2311 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2316 TString triggerClass(trigger->String());
2320 while ( ( str = static_cast<TObjString*>(nextFile())) )
2322 TObjArray* a = static_cast<TObjArray*>(m->GetValue(str->String().Data()));
2324 AliAnalysisMuMuJpsiResult* r = static_cast<AliAnalysisMuMuJpsiResult*>(a->FindObject(triggerClass.Data()));
2328 int aodtrain,esdpass,runnumber;
2330 if ( DecodeFileName(str->String().Data(),period,esdpass,aodtrain,runnumber) )
2332 runnumbers.insert(runnumber);
2334 runMin = TMath::Min(runMin,runnumber);
2335 runMax = TMath::Max(runMax,runnumber);
2337 x_jpsirate[triggerClass.Data()].push_back(runnumber);
2338 xerr_jpsirate[triggerClass.Data()].push_back(0.5);
2342 TString what("RateJpsi");
2345 what = "AccEffJpsi";
2348 if ( TMath::Finite(r->GetValue("SigmaJpsi")) && r->NofTriggers() > 10 )
2350 y = 100*r->GetValue(what.Data());
2351 yerr = 100*r->GetErrorStat(what.Data());
2353 if (!efficiencies.empty() )
2355 if (efficiencies.count(runnumber))
2357 y /= ( efficiencies[runnumber] );
2365 ymin = TMath::Min(ymin,y);
2366 ymax = TMath::Max(ymax,y);
2368 sumw += y*r->NofTriggers();
2369 n += r->NofTriggers();
2372 y_jpsirate[triggerClass.Data()].push_back(y);
2373 yerr_jpsirate[triggerClass.Data()].push_back(yerr);
2377 AliInfoClass(Form("Trigger %30s ponderated mean is %7.2f",trigger->String().Data(),sumw/n));
2382 TString canvasName("cJpsiRateEvolution");
2384 if ( !efficiencies.empty() )
2386 canvasName += "Corr";
2389 TCanvas* c = new TCanvas(canvasName.Data(),canvasName.Data());
2393 Int_t nbins = runnumbers.size();
2397 if ( CompactGraphs() )
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);
2406 gStyle->SetOptTitle(0);
2407 gStyle->SetOptStat(0);
2409 if (!CompactGraphs())
2411 h->GetXaxis()->SetNoExponent();
2415 std::set<int>::const_iterator it;
2418 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2420 h->GetXaxis()->SetBinLabel(i,Form("%d",*it));
2423 h->GetXaxis()->SetNdivisions(1,kFALSE);
2431 AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath);
2432 ts.DrawFills(ymin,ymax);
2435 h->Draw("sameaxis");
2437 //c->RedrawAxis("g");
2439 nextTrigger.Reset();
2442 int color[] = { 2,1,4,5,6 };
2443 int marker[] = { 20,23,25,21,22 };
2445 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
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()];
2452 TGraphErrors* g = new TGraphErrors(x.size(),&x[0],&y[0],&xerr[0],&yerr[0]);
2454 g->SetLineColor(1);//color[i]);
2455 g->SetMarkerColor(color[i]);
2456 g->SetMarkerStyle(marker[i]);
2457 g->SetMarkerSize(0.7);
2458 g->GetXaxis()->SetNoExponent();
2460 if ( CompactGraphs() )
2462 AliAnalysisMuMuGraphUtil::Compact(*g);
2466 TString gname(trigger->String());
2467 gname.ReplaceAll("-","_");
2468 g->SetName(gname.Data());
2471 Double_t m2 = g->GetMean(2);
2473 TLine* line = new TLine(runMin,m2,runMax,m2);
2474 line->SetLineColor(color[i]);
2477 AliInfoClass(Form("TRIGGER %s MEAN %7.2f",trigger->String().Data(),m2));
2484 //_____________________________________________________________________________
2485 TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const char* event1,
2486 const char* trigger2, const char* event2,
2488 Bool_t asRejection) const
2490 if (!CC()) return 0x0;
2492 const std::set<int>& runnumbers = RunNumbers();
2494 TGraphErrors* g = new TGraphErrors(runnumbers.size());
2496 std::set<int>::const_iterator it;
2499 Double_t ymin(TMath::Limits<double>::Max());
2500 Double_t ymax(TMath::Limits<double>::Min());
2502 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
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));
2511 if ( fCorrectionPerRun )
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
2516 // FIXME: should get the correction error here
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);
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);
2537 gStyle->SetOptStat(0);
2541 hframe->GetXaxis()->SetNoExponent();
2543 hframe->GetYaxis()->SetTitle(asRejection ? "Rejection (%)" : "Ratio");
2546 g->SetTitle(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2));
2547 g->GetXaxis()->SetNoExponent();
2550 AliAnalysisTriggerScalers ts(RunNumbers(),fgOCDBPath.Data());
2554 ts.DrawFills(ymin,ymax);
2559 std::map<std::string, std::pair<int,int> > periods;
2561 ts.GetLHCPeriodBoundaries(periods);
2563 TLegend* legend = new TLegend(0.15,0.82,0.90,0.92);
2564 legend->SetFillColor(0);
2568 for ( std::map<std::string, std::pair<int,int> >::const_iterator pit = periods.begin(); pit != periods.end(); ++pit )
2570 std::string period = pit->first;
2571 int run1 = (pit->second).first;
2572 int run2 = (pit->second).second;
2574 for ( std::set<int>::const_iterator rit = RunNumbers().begin(); rit != RunNumbers().end(); ++ rit )
2576 if ( (*rit) >= run1 && (*rit) <= run2 )
2581 AliInfo(Form("Period %s runs %6d-%6d ; %d actual runs",period.c_str(),run1,run2,nruns));
2583 g->Fit("pol0","+Q","",run1,run2);
2584 TF1* func = static_cast<TF1*>(g->GetListOfFunctions()->Last());
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)));
2594 legend->SetNColumns(3);
2596 Double_t mean = TMath::Mean(g->GetN(),g->GetY());
2597 Double_t rms = TMath::RMS(g->GetN(),g->GetY());
2599 legend->AddEntry("",Form("Mean %5.2f RMS %5.2f (%5.2f %%)",mean,rms,(mean) ? 100.0*rms/mean : 0.0),"");
2606 //_____________________________________________________________________________
2607 void AliAnalysisMuMu::Print(Option_t* opt) const
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;
2617 std::cout << RunNumbers().size() << " runs";
2618 if ( fCorrectionPerRun )
2620 std::cout << " with correction factors";
2622 std::cout << std::endl;
2624 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2627 if ( fCorrectionPerRun )
2629 std::cout << Form("(%e)",fCorrectionPerRun->GetY()[i]);
2634 std::cout << std::endl;
2639 if ( sopt.Contains("BIN") && BIN() )
2641 std::cout << "Binning : " << std::endl;
2643 topt.ReplaceAll("BIN","");
2644 BIN()->Print(topt.Data());
2646 if ( sopt.Contains("MC") && MC() )
2649 topt.ReplaceAll("MC","");
2650 MC()->Print(topt.Data());
2652 if ( sopt.Contains("CC") && CC() )
2654 CC()->Print("trigger/event");
2657 if ( sopt.Contains("SIZE") )
2659 TFile* f = ReOpen(fFilename,"READ");
2660 TIter next(f->GetListOfKeys());
2663 while ( ( key = static_cast<TKey*>(next()) ) )
2665 std::cout << key->GetName() << " " << key->GetNbytes() << " " << key->GetObjlen() << std::endl;
2670 //_____________________________________________________________________________
2672 AliAnalysisMuMu::ReadFileList(const char* filelist)
2675 // read the filelist and try to order it by runnumber (or periods)
2677 // filelist can either be a real filelist (i.e. a text file containing
2678 // root filenames) or a root file itself.
2683 TObjArray* files = new TObjArray;
2684 files->SetOwner(kTRUE);
2686 TString sfilelist(ExpandPathName(filelist));
2688 if ( sfilelist.EndsWith(".root") )
2690 files->Add(new TObjString(sfilelist.Data()));
2694 std::set<std::string> sorting;
2695 std::map<std::string,std::string> filemap;
2697 std::ifstream in(sfilelist.Data());
2700 int aodtrain,esdpass,runnumber;
2702 while ( in.getline(line,1022,'\n') )
2704 TString sline(line);
2705 if (sline.BeginsWith("#")) continue;
2707 DecodeFileName(line,period,esdpass,aodtrain,runnumber);
2709 AliDebugClass(1,Form("line %s => period %s esdpass %d aodtrain %d runnumber %09d",
2710 line,period.Data(),esdpass,aodtrain,runnumber));
2712 TString key(Form("%d",runnumber));
2714 if ( runnumber <= 0 )
2718 sorting.insert(key.Data());
2719 filemap.insert(std::make_pair<std::string,std::string>(key.Data(),line));
2724 std::set<std::string>::const_iterator it;
2726 for ( it = sorting.begin(); it != sorting.end(); ++it )
2728 files->Add(new TObjString(filemap[*it].c_str()));
2733 //_____________________________________________________________________________
2734 TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
2736 /// Tries to reopen the file with a new mode
2738 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(filename));
2745 f = TFile::Open(filename,mode);
2747 if ( !f || !f->IsOpen() )
2749 AliError(Form("Cannot open file %s in mode %s",filename,mode));
2756 //_____________________________________________________________________________
2758 AliAnalysisMuMu::RunNumberFromFileName(const char* filename)
2761 int esdpass,aodtrain,runnumber;
2762 Bool_t ok = DecodeFileName(filename,period,esdpass,aodtrain,runnumber);
2763 if ( ok ) return runnumber;
2767 //_____________________________________________________________________________
2768 void AliAnalysisMuMu::SetColorScheme()
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");
2774 gStyle->SetGridColor(AliAnalysisMuMu::kBlue);
2776 gStyle->SetFrameLineColor(AliAnalysisMuMu::kBlue);
2777 gStyle->SetAxisColor(AliAnalysisMuMu::kBlue,"xyz");
2778 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue,"xyz");
2780 gStyle->SetTitleColor(AliAnalysisMuMu::kBlue);
2781 gStyle->SetTitleTextColor(AliAnalysisMuMu::kBlue);
2782 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue);
2783 gStyle->SetStatTextColor(AliAnalysisMuMu::kBlue);
2785 gStyle->SetOptStat(0);
2788 //_____________________________________________________________________________
2789 Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
2791 /// Sets the graph used to correct values per run
2792 delete fCorrectionPerRun;
2793 fCorrectionPerRun=0x0;
2795 // check that corr has the same runs as we do
2799 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2801 Int_t corrRun = TMath::Nint(corr.GetX()[i]);
2804 AliError(Form("%d-th run mistmatch %d vs %d",i,corrRun,*it));
2811 fCorrectionPerRun = new TGraphErrors(corr.GetN());
2813 TFormula* tformula(0x0);
2814 if ( strlen(formula) > 0 )
2816 tformula = new TFormula("SetCorrectionPerRunFormula",formula);
2821 for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2823 Double_t y = corr.GetY()[i];
2827 y = tformula->Eval(y);
2829 fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
2838 //_____________________________________________________________________________
2839 void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuJpsiResult& r)
2841 /// Set the "NofInput" variable(s) of one result
2843 TString hname(Form("MinvUS%s",r.Bin().AsString().Data()));
2845 TH1* hinput = fMergeableCollection->Histo("/INPUT/INYRANGE",hname.Data());
2849 AliError(Form("Got a simulation file where I did not find histogram /INPUT/INYRANGE/%s",hname.Data()));
2854 r.SetNofInputParticles(*hinput);
2858 //_____________________________________________________________________________
2859 AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
2861 /// Shortcut method to get to a spectra
2862 if (!MC()) return 0x0;
2864 return static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(fullpath));
2867 //_____________________________________________________________________________
2868 void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
2870 Bool_t orderByTriggerCount)
2872 // Give the fraction of triggers (in triggerList) relative
2873 // to what is expected in the scalers
2875 TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
2877 AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
2879 AliLog::SetGlobalLogLevel(AliLog::kFatal);
2881 if (!fMergeableCollection || !fCounterCollection) return;
2883 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2884 TIter nextRun(runs);
2886 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2887 TIter nextTrigger(triggers);
2890 TObjString* strigger;
2892 TString striggerList(triggerList);
2895 ULong64_t totalExpected(0);
2897 std::multimap<ULong64_t,std::string> messages;
2899 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2901 msg.Form("RUN %09d ",srun->String().Atoi());
2910 nextTrigger.Reset();
2912 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2914 if ( !striggerList.Contains(strigger->String().Data()) )
2918 ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2919 strigger->String().Data(),"ALL",srun->String().Atoi())));
2921 ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
2924 nmax = TMath::Max(n,nmax);
2927 totalExpected += expected;
2929 msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
2930 (n>expected ? "!" : " "));
2932 if ( expected > 0 ) {
2933 msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
2943 if (!orderByTriggerCount)
2945 std::cout << msg.Data() << std::endl;
2949 messages.insert(std::make_pair(nmax,static_cast<std::string>(msg.Data())));
2954 std::multimap<ULong64_t,std::string>::const_reverse_iterator it;
2956 ULong64_t current(0);
2959 for ( it = messages.rbegin(); it != messages.rend(); ++it )
2962 current += it->first;
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;
2967 std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
2968 total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
2971 AliLog::SetGlobalLogLevel(oldLevel);
2976 //_____________________________________________________________________________
2977 void AliAnalysisMuMu::UnsetCorrectionPerRun()
2979 // drop the correction factors
2980 delete fCorrectionPerRun;
2981 fCorrectionPerRun=0x0;
2984 //_____________________________________________________________________________
2985 void AliAnalysisMuMu::Update()
2987 /// update the current file with memory
2989 if (!CC() || !MC()) return;
2991 ReOpen(fFilename,"UPDATE");
2995 MC()->Write("MC",TObject::kSingleKey|TObject::kOverwrite);
2998 ReOpen(fFilename,"READ");
3000 GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
3003 //_____________________________________________________________________________
3004 Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
3007 AliAnalysisMuMu m(filename);
3012 //_____________________________________________________________________________
3013 Bool_t AliAnalysisMuMu::Upgrade()
3015 /// Upgrade the current file
3016 /// - from single list to one key per object, if needed
3017 /// - from histogramCollection to mergeableCollection, if needed
3019 TFile* f = ReOpen(fFilename,"UPDATE");
3021 TList* list = static_cast<TList*>(f->Get("chist"));
3025 // really old file where everything was in a single list
3027 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(list->At(0));
3028 AliCounterCollection* cc = static_cast<AliCounterCollection*>(list->At(1));
3030 AliMergeableCollection* mc = hc->Convert();
3034 mc->Write("MC",TObject::kSingleKey);
3035 cc->Write("CC",TObject::kSingleKey);
3037 f->Delete("chist;*");
3044 AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(f->Get("HC"));
3048 // old file with histogram collection instead of mergeable collection
3050 AliMergeableCollection* mc = hc->Convert();
3054 mc->Write("MC",TObject::kSingleKey);
3067 //_____________________________________________________________________________
3068 void AliAnalysisMuMu::ComputeFnorm()
3070 /// Compute the CMUL to CINT ratio(s)
3074 MC()->Prune("/FNORM");
3076 AliAnalysisMuMuFnorm computer(*(CC()),AliAnalysisMuMuFnorm::kMUL,fgOCDBPath.Data(),fgIsCompactGraphs);
3078 computer.ComputeFnorm();
3080 AliMergeableCollection* fnorm = computer.DetachMC();
3082 MC()->Attach(fnorm,"/FNORM/");
3087 //_____________________________________________________________________________
3088 AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
3090 /// Correct one spectra
3094 AliError("Cannot compute corrected yield without associated MC file !");
3098 const char* accEffSubResultName="COUNTJPSI:1";
3100 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
3101 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
3105 AliError("could not get real spectra");
3111 AliError("could not get sim spectra");
3115 realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
3122 //_____________________________________________________________________________
3123 AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
3127 AliError("Cannot compute corrected yield without associated MC file !");
3131 const char* accEffSubResultName="COUNTJPSI:1";
3133 AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
3137 AliError("could not get real spectra");
3141 if (!realSpectra->HasValue("CoffNofJpsi"))
3143 if (!CorrectSpectra(type,flavour))
3145 AliError("Could not get corrected spectra");
3150 AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
3154 AliErrorClass("could not get sim spectra");
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"));
3162 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
3163 TObjArray* bins = binning->CreateBinObjArray();
3164 TIter nextBin(bins);
3165 AliAnalysisMuMuBinning::Range* bin;
3167 AliAnalysisMuMuJpsiResult* r;
3169 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
3171 r = static_cast<AliAnalysisMuMuJpsiResult*>(realSpectra->BinContentArray()->At(i));
3173 StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
3175 AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
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));
3182 r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3183 nofCINT7,TMath::Sqrt(nofCINT7)));
3185 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
3187 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
3190 r->Set("YJpsi",yield,yieldError);
3192 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
3193 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
3205 //_____________________________________________________________________________
3206 AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* realFile, const char* simFile, const char* type,
3207 const char* direction)
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
3215 const char* accEffSubResultName="COUNTJPSI:1";
3217 TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
3218 ydist.SetParameter(0,1.);
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)));
3228 AliAnalysisMuMu real(realFile);
3229 AliAnalysisMuMu sim(simFile);
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)));
3237 AliErrorClass("could not get real spectra");
3243 AliErrorClass("could not get sim spectra");
3247 AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
3248 corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
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"));
3254 AliAnalysisMuMuBinning* binning = realSpectra->Binning();
3255 TObjArray* bins = binning->CreateBinObjArray();
3256 TIter nextBin(bins);
3257 AliAnalysisMuMuBinning::Range* bin;
3259 AliAnalysisMuMuJpsiResult* r;
3261 Int_t n = bins->GetLast();
3263 TObjArray finalBins(n+1);
3264 finalBins.SetOwner(kTRUE);
3266 TObjArray finalResults(n+1);
3267 finalResults.SetOwner(kFALSE);
3269 while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
3271 Double_t ylowlab = bin->Xmin();
3272 Double_t yhighlab = bin->Xmax();
3274 Double_t ylowcms, yhighcms;
3275 Double_t ylownorm, yhighnorm;
3277 if ( bin->IsNullObject() )
3283 if ( strcmp(direction,"pPb")==0 )
3285 ylowcms = TMath::Abs(yhighlab) - rapidityShift;
3286 yhighcms = TMath::Abs(ylowlab) - rapidityShift;
3287 ylownorm = ylowcms/ymax;
3288 yhighnorm = yhighcms/ymax;
3292 ylowcms = ylowlab - rapidityShift;
3293 yhighcms = yhighlab - rapidityShift;
3294 ylownorm = -yhighcms/ymax;
3295 yhighnorm = -ylowcms/ymax;
3299 Double_t brsigmapp = ydist.Integral(ylownorm,yhighnorm);
3300 Double_t brsigmappError = 0.0; // FIXME
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));
3305 r = static_cast<AliAnalysisMuMuJpsiResult*>(corrSpectra->BinContentArray()->At(i)->Clone());
3307 AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
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));
3314 r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3315 nofCINT7,TMath::Sqrt(nofCINT7)));
3317 Double_t yield = r->GetValue("CorrNofJpsi") * mbeq;
3319 Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
3322 r->Set(Form("Y%sJpsi",direction),yield,yieldError);
3324 Double_t raa = yield/(tab*brsigmapp);
3325 Double_t raaError = AliAnalysisMuMuResult::ErrorABC(yield,yieldError,
3327 brsigmapp,brsigmappError);
3328 r->Set(Form("R%sJpsi",direction),raa,raaError);
3330 r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
3331 r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
3333 AliAnalysisMuMuBinning::Range* bincm = new AliAnalysisMuMuBinning::Range(bin->Particle(),bin->Type(),ylowcms,yhighcms);
3337 finalBins.Add(bincm);
3338 finalResults.Add(r);
3345 AliAnalysisMuMuSpectra* spectra = new AliAnalysisMuMuSpectra(type,direction);
3347 for ( i = 0; i <= n; ++i )
3350 if ( strcmp(direction,"pPb")==0 )
3355 r = static_cast<AliAnalysisMuMuJpsiResult*>(finalResults.At(j));
3357 bin = static_cast<AliAnalysisMuMuBinning::Range*>(finalBins.At(j));
3359 spectra->AdoptResult(*bin,r);
3368 //_____________________________________________________________________________
3369 void AliAnalysisMuMu::ComputeJpsi(const char* runlist, const char* prefix, const char* what, const char* binningFlavour)
3371 /// Call the Jpsi method for each file
3373 std::vector<int> runs;
3374 AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3376 for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3378 Int_t runNumber = runs[i];
3380 TString filename(Form("RUNBYRUN/%s_%09d.saf.root",prefix,runNumber));
3382 if ( gSystem->AccessPathName(filename.Data())==kTRUE ) continue;
3384 AliAnalysisMuMu m(filename.Data());
3385 m.Jpsi(what,binningFlavour);
3389 //_____________________________________________________________________________
3390 AliAnalysisMuMuSpectra*
3391 AliAnalysisMuMu::CombineSpectra(const char* runlist,
3392 const char* realPrefix,
3393 const char* simPrefix,
3394 const char* quantity,
3395 const char* variable,
3396 const char* binningFlavour)
3398 std::vector<int> runs;
3399 AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3401 std::vector<double> vw;
3402 std::vector<AliAnalysisMuMuSpectra*> vspectra;
3404 for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3406 Int_t runNumber = runs[i];
3408 TString filename(Form("RUNBYRUN/%s_%09d.saf.root",realPrefix,runNumber));
3410 AliAnalysisMuMu mreal(filename.Data());
3412 std::cout << filename.Data() << std::flush << std::endl;
3414 if ( !mreal.CC() || !mreal.MC() )
3416 AliErrorClass(Form("Have to skip %s CC=%p MC=%p",filename.Data(),mreal.CC(),mreal.MC()));
3420 TGraph* g = static_cast<TGraph*>(mreal.MC()->GetObject("/FNORM/GRAPHS/NMBeqBest2"));
3424 mreal.ComputeFnorm();
3425 g = static_cast<TGraph*>(mreal.MC()->GetObject("/FNORM/GRAPHS/NMBeqBest2"));
3430 if ( TMath::Nint(g->GetX()[0]) != runNumber )
3432 AliErrorClass("Wrong run number in NMBeqBest2 graph !");
3436 filename.Form("RUNBYRUN/%s_%09d.saf.root",simPrefix,runNumber);
3438 AliAnalysisMuMu msim(filename.Data());
3440 TString spectraName(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",variable));
3441 if (strlen(binningFlavour))
3444 spectraName += binningFlavour;
3447 AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject(spectraName.Data()));
3451 msim.Jpsi(quantity,binningFlavour);
3453 s = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject(spectraName.Data()));
3458 AliErrorClass(Form("Could not find spectra for sim %s",filename.Data()));
3462 vw.push_back(g->GetY()[0]);
3463 vspectra.push_back(static_cast<AliAnalysisMuMuSpectra*>(s->Clone()));
3468 for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3471 vspectra[i]->Print();
3474 // normalize the weights
3476 list.SetOwner(kTRUE);
3478 for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3482 vspectra[i]->SetWeight(vw[i]);
3484 if ( i > 0 ) list.Add(vspectra[i]);
3487 std::cout << "before merging" << std::endl;
3488 for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3490 std::cout << " ---------------- " << i << std::endl;
3492 vspectra[i]->Print();
3495 vspectra[0]->Merge(&list);
3502 //_____________________________________________________________________________
3503 TGraph* AliAnalysisMuMu::ResultEvolution(const char* runlist, const char* realPrefix, const char* simPrefix, const char* what, Bool_t forceRecomputation)
3505 std::vector<int> runs;
3506 AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3508 TGraphErrors* g = new TGraphErrors(runs.size());
3510 TString direction("Pbp");
3511 TString check(realPrefix);
3514 if (check.Contains("LHC13B") ||
3515 check.Contains("LHC13C") ||
3516 check.Contains("LHC13D") ||
3517 check.Contains("LHC13E")
3523 AliInfoClass(Form("direction %s",direction.Data()));
3528 TString subResultName("");
3530 TString swhat(what);
3532 TObjArray* parts = swhat.Tokenize(":");
3534 if ( parts->GetEntries() > 1 )
3536 subResultName = static_cast<TObjString*>(parts->At(1))->String();
3541 for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3543 Int_t runNumber = runs[i];
3545 AliDebugClass(1,Form("RUN %09d",runNumber));
3547 TString realFile(Form("RUNBYRUN/%s_%09d.saf.root",realPrefix,runNumber));
3549 TString simFileName(Form("RUNBYRUN/%s_%09d.saf.root",simPrefix,runNumber));
3551 TString simFile(simFileName);
3553 TString resultName(Form("%s%sJpsi",what,direction.Data()));
3555 if ( swhat == "Fnorm")
3557 resultName = "Fnorm";
3559 else if ( swhat.Contains("Acc") )
3561 resultName.ReplaceAll(direction,"");
3564 AliAnalysisMuMu* mreal(0x0);
3566 if ( !swhat.Contains("Acc"))
3568 mreal = new AliAnalysisMuMu(realFile);
3571 AliAnalysisMuMuSpectra* real(0x0);
3575 real = static_cast<AliAnalysisMuMuSpectra*>(mreal->MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3577 if (!real || forceRecomputation)
3580 real = static_cast<AliAnalysisMuMuSpectra*>(mreal->MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3583 AliErrorClass(Form("Could not get real spectra for run %d",runNumber));
3590 AliAnalysisMuMu msim(simFile);
3592 AliAnalysisMuMuSpectra* sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3594 if (!sim || forceRecomputation)
3596 msim.SetEventSelectionList("ALL");
3598 sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3601 AliErrorClass(Form("Could not get sim spectra for run %d",runNumber));
3607 AliAnalysisMuMuSpectra* corrected = 0x0;
3611 corrected = AliAnalysisMuMu::RABy(realFile.Data(),simFile.Data(),"",direction.Data());
3618 AliAnalysisMuMuResult* result = static_cast<AliAnalysisMuMuResult*>(corrected->BinContentArray()->First());
3622 Double_t value = result->GetValue(resultName.Data(),subResultName.Data());
3623 Double_t error = result->GetErrorStat(resultName.Data(),subResultName.Data());
3625 g->SetPoint(i,runNumber,value);
3626 g->SetPointError(i,1,error);
3632 n = mreal->CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL/run:%d",runNumber));
3636 if ( error > 0 ) w = 1.0/error/error;
3641 // std::cout << result->SubResults() << " " << result->GetError(resultName.Data()) << std::endl;
3646 gStyle->SetOptFit(1);
3650 g->GetXaxis()->SetTitle("Run number");
3651 g->GetXaxis()->SetNoExponent();
3652 if ( TString(what) == "Y" )
3654 g->GetYaxis()->SetTitle("J/#psi yield");
3656 else if ( TString(what) == "R" )
3658 g->GetYaxis()->SetTitle(Form("R_{%s}^{J/#psi}",direction.Data()));
3660 else if ( TString(what).Contains("Acc") )
3662 g->GetYaxis()->SetTitle("Acc#timesEff_{J/#psi}");
3664 else if ( TString(what).Contains("Fnorm") )
3666 g->GetYaxis()->SetTitle("CINT7 / CINT7&0MUL");
3669 if (CompactGraphs())
3671 AliAnalysisMuMuGraphUtil::Compact(*g);
3676 AliInfoClass(Form("Weighted Mean %e",mean));