]> 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
Always delete TObjArrays created by TString::Tokenize (Ruben)
[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
21#include "AliAnalysisTriggerScalers.h"
22#include "AliCounterCollection.h"
23#include "AliHistogramCollection.h"
24#include "AliLog.h"
25#include "Riostream.h"
26#include "TArrayL64.h"
27#include "TAxis.h"
28#include "TCanvas.h"
29#include "TF1.h"
30#include "TFile.h"
31#include "TGraphErrors.h"
32#include "TGrid.h"
33#include "TH1.h"
34#include "TH2.h"
35#include "TLegend.h"
36#include "TLegendEntry.h"
37#include "TList.h"
38#include "TMath.h"
39#include "TObjArray.h"
40#include "TObjString.h"
41#include "TROOT.h"
42#include "TStyle.h"
43#include "TSystem.h"
44#include <map>
45#include <set>
46#include <string>
47#include "TParameter.h"
48#include "TMap.h"
49#include "TFitResult.h"
50#include "TLine.h"
51#include "TASImage.h"
52#include "TPaveText.h"
53#include "TStyle.h"
54#include "TColor.h"
55
56ClassImp(AliAnalysisMuMu)
57
58void Add(TObjArray* a, AliAnalysisMuMu::Result* r)
59{
60 if ( r ) a->Add(r);
61}
62
63void ALICEseal(Double_t xPad, Double_t yPad, TList* extraLines)
64{
65 TVirtualPad* currPad = gPad;
66 TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.76,0.7,0.90,0.87);
67 myPadLogo->SetFillColor(kWhite);
68 myPadLogo->SetFillStyle(1001);
69 myPadLogo->SetBorderMode(0);
70 myPadLogo->SetBorderSize(2);
71 myPadLogo->SetFrameBorderMode(0);
72 myPadLogo->SetLeftMargin(0.0);
73 myPadLogo->SetTopMargin(0.0);
74 myPadLogo->SetBottomMargin(0.0);
75 myPadLogo->SetRightMargin(0.0);
76 myPadLogo->Draw();
77 myPadLogo->cd();
78 TASImage *myAliceLogo = new TASImage("$HOME/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.gif");
79 myAliceLogo->Draw();
80 currPad->cd();
81 Double_t x1 = xPad - 0.07, y1 = yPad - 0.06;
82 Double_t x2 = x1 + 0.25, y2 = y1 + 0.08;
83 TPaveText* t2=new TPaveText(x1+0.06,y1-0.06,x2-0.06,y2-0.06,"NDC");
84 // t2->SetFillStyle(0);
85 t2->SetFillColor(kWhite);
86 t2->SetBorderSize(0);
87 t2->SetTextColor(kBlack);
88 t2->SetTextFont(52);
89 t2->SetTextSize(0.035);
90 TDatime dt;
91 TString today = Form("%02i/%02i/%4i", dt.GetDay(), dt.GetMonth(), dt.GetYear());
92 t2->AddText(0.,0.,today.Data());
93 t2->Draw();
94
95 if ( extraLines )
96 {
97 int n = extraLines->GetSize();
98 TPaveText* t3 = new TPaveText(xPad-0.07,yPad-0.06*2,xPad-0.07+0.25,yPad-0.06*(n+3),"NDC");
99 t3->SetFillColor(kWhite);
100 t3->SetBorderSize(0);
101 t3->SetTextColor(kBlack);
102 t3->SetTextFont(42);
103 t3->SetTextSize(0.035);
104
105 TIter next(extraLines);
106 TObjString* str;
107
108 while ( ( str = static_cast<TObjString*>(next()) ) )
109 {
110 t3->AddText(str->String().Data());
111 }
112 t3->SetFillColor(kWhite);
113 t3->Draw();
114 }
115}
116
117Double_t funcCB(Double_t* xx, Double_t* par)
118{
119 Double_t N = par[0];
120 Double_t alpha = par[1];
121 Double_t n = par[2];
122 Double_t mean = par[3];
123 Double_t sigma = par[4];
124
125 Double_t x = xx[0];
126
127 Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
128 Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
129
130 Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
131
132 if ( y > alpha*-1.0 )
133 {
134 return N*TMath::Exp(-0.5*y*y);
135 }
136 else
137 {
138 return N*A*TMath::Power(B-y,-n);
139 }
140}
141
142Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
143{
144 Double_t x = xx[0];
145
146 Double_t g = par[0]*TMath::Gaus(x,par[1],par[2]);
147
148 Double_t jpsi = funcCB(xx,par+3);
149
150 Double_t expo = par[8]*TMath::Exp(par[9]*x);
151
152 return g+expo+jpsi;
153}
154
155Double_t funcJpsiJpsiPrimeCustom(Double_t* xx, Double_t* par)
156{
157 Double_t N = par[0];
158 Double_t alpha = par[1];
159 Double_t n = par[2];
160 Double_t mean = par[3];
161 Double_t sigma = par[4];
162 Double_t alphaprime = par[5];
163 Double_t nprime = par[6];
164
165 Double_t x = xx[0];
166
167 Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
168 Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
169 Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
170 Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
171
172 Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
173
174 Double_t cb(0);
175
176 if ( y > alphaprime )
177 {
178 cb = N*C*TMath::Power(D+y,-nprime);
179 }
180 else if ( y > alpha*-1.0 )
181 {
182 cb = N*TMath::Exp(-0.5*y*y);
183 }
184 else
185 {
186 cb = N*A*TMath::Power(B-y,-n);
187 }
188
189 if ( x < mean )
190 {
191 return cb + par[7] + par[8]*x; // gaus + pol1
192 }
193 else
194 {
195 Double_t yprime = (x-par[10])/par[11];
196 return cb + par[9]*TMath::Exp(-0.5*yprime*yprime)+par[12]*TMath::Exp(-par[13]*x);
197 // gaus (j/psi) + gaus (psi') + expo
198 }
199}
200
201
202Double_t funcCB2(Double_t* xx, Double_t* par)
203{
204 Double_t N = par[0];
205 Double_t alpha = par[1];
206 Double_t n = par[2];
207 Double_t mean = par[3];
208 Double_t sigma = par[4];
209 Double_t alphaprime = par[5];
210 Double_t nprime = par[6];
211
212 Double_t x = xx[0];
213
214 Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
215 Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
216 Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
217 Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
218
219 Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
220
221 if ( y > alphaprime )
222 {
223 return N*C*TMath::Power(D+y,-nprime);
224 }
225 else if ( y > alpha*-1.0 )
226 {
227 return N*TMath::Exp(-0.5*y*y);
228 }
229 else
230 {
231 return N*A*TMath::Power(B-y,-n);
232 }
233}
234
235Double_t funcJpsiJpsiPrime(Double_t* xx, Double_t* par)
236{
237 Double_t jpsi = funcCB(xx,par);
238 Double_t jpsiprime = funcCB2(xx,par+5);
239
240 int n = 10;
241 Double_t x = xx[0];
242
243 Double_t e1 = par[n]*TMath::Exp(par[n+1]*x);
244 Double_t e2 = par[n+2]*TMath::Exp(par[n+3]*x);
245
246 Double_t e = e1;
247
248 if ( x > par[3] ) e=e2;
249
250 return jpsi+jpsiprime+e;
251}
252
253Double_t funcJpsiCBE(Double_t* xx, Double_t* par)
254{
255 // CB + expo
256
257 Double_t jpsi = funcCB(xx,par);
258
259 Double_t x = xx[0];
260
261 Double_t e1 = par[5]*TMath::Exp(par[6]*x);
262
263 return jpsi+e1;
264}
265
266
267Double_t funcJpsiPCBE(Double_t* xx, Double_t* par)
268{
269 Double_t x = xx[0];
270
271 Double_t pol2 = par[0] + par[1]*x + par[2]*x*x;
272
273 Double_t jpsi = funcCB(xx,par+3);
274
275 Double_t expo = par[8]*TMath::Exp(par[9]*x);
276
277 return pol2+jpsi+expo;
278}
279
280Double_t funcJpsiECBE(Double_t* xx, Double_t* par)
281{
282 // CB + expo
283
284 Double_t jpsi = funcCB(xx,par+2);
285
286 Double_t x = xx[0];
287
288 Double_t e1 = par[0]*TMath::Exp(par[1]*x);
289
290 Double_t e2 = par[7]*TMath::Exp(par[8]*x);
291
292 return e1+e2+jpsi;
293}
294
295const char* NormalizeName(const char* name, const char* suffix)
296{
297 TString str(Form("%s_%s",name,suffix));
298
299 str.ReplaceAll("-","_");
300 str.ReplaceAll("/","%");
301
302 return str.Data();
303}
304
305//_____________________________________________________________________________
306AliAnalysisMuMu::Result::~Result()
307{
308 delete fHC;
309 delete fMap;
310}
311
312//_____________________________________________________________________________
313void AliAnalysisMuMu::Result::FitJpsiJpsiPrimeCustom(TH1& h)
314{
315 std::cout << "Fit with jpsi + jpsiprime' (custom)" << std::endl;
316
317 const Double_t xmin(1.5);
318 const Double_t xmax(8.0);
319
320 fFitTotal = new TF1("fFitTotal",funcJpsiJpsiPrimeCustom,xmin,xmax,14);
321 fFitTotal->SetLineColor(4);
322
323 fFitTotal->SetParName(0,"cstecb");
324 fFitTotal->SetParName(1,"alpha");
325 fFitTotal->SetParName(2,"n");
326 fFitTotal->SetParName(3,"meanjpsi");
327 fFitTotal->SetParName(4,"sigmajpsi");
328 fFitTotal->SetParName(5,"alphaprime");
329 fFitTotal->SetParName(6,"nprime");
330 fFitTotal->SetParName(7,"cstepol1");
331 fFitTotal->SetParName(8,"slopepol1");
332 fFitTotal->SetParName(9,"cstegaus");
333 fFitTotal->SetParName(10,"meanpsiprime");
334 fFitTotal->SetParName(11,"sigmapsiprime");
335 fFitTotal->SetParName(12,"csteexpo");
336 fFitTotal->SetParName(13,"slopeexpo");
337
338 fFitTotal->SetParameter( 0,1);
339
340 const char* fitOption = "SQBR+";
341 const Double_t alphaMC = 0.936;
342 const Double_t nMC = 4.44;
343 const Double_t alphaprimeMC = 1.60;
344 const Double_t nprimeMC = 3.23;
345
346 TF1* fcb = new TF1("cb",funcCB2,2.9,3.3,7);
347 fcb->SetParameters(1,1.0,4.0,3.1,0.1,1.5,3);
348
349 fcb->SetParLimits(3,3,4);
350 fcb->SetParLimits(4,0,1);
351
352 fcb->FixParameter(1,alphaMC);
353 fcb->FixParameter(2,nMC);
354 fcb->FixParameter(5,alphaprimeMC);
355 fcb->FixParameter(6,nprimeMC);
356
357 TFitResultPtr rcb = h.Fit(fcb,fitOption,"",2.9,3.3);
358
359 if (!rcb.Get())
360 {
361 return;
362 }
363
364 fFitTotal->SetParameter(0,rcb->Parameter(0));
365 fFitTotal->SetParameter(1,rcb->Parameter(1)); fFitTotal->SetParLimits(1,0,10); // alpha
366 fFitTotal->SetParameter(2,rcb->Parameter(2)); fFitTotal->SetParLimits(2,1,10); // n
367 fFitTotal->SetParameter(3,rcb->Parameter(3)); fFitTotal->SetParLimits(3,3.0,3.5); // mean
368 fFitTotal->SetParameter(4,rcb->Parameter(4)); fFitTotal->SetParLimits(4,0,1); // sigma
369 fFitTotal->SetParameter(5,rcb->Parameter(5)); fFitTotal->SetParLimits(1,0,10); // alphaprime
370 fFitTotal->SetParameter(6,rcb->Parameter(6)); fFitTotal->SetParLimits(2,1,10); // nprime
371
372 fFitTotal->FixParameter(1,alphaMC);
373 fFitTotal->FixParameter(2,nMC);
374 fFitTotal->FixParameter(5,alphaprimeMC);
375 fFitTotal->FixParameter(6,nprimeMC);
376
377 TF1* fge = new TF1("fge","gaus(0)+expo(3)",3.5,4.4);
378 fge->SetParameters(1,3.6,0.25,1,1);
379 TFitResultPtr rpsiprime = h.Fit(fge,fitOption,"",3.5,4.4);
380
381 if (static_cast<int>(rpsiprime))
382 {
383 AliInfo("Will fix psiprime parameters");
384 fFitTotal->FixParameter(9,0);
385 fFitTotal->FixParameter(10,3.7);
386 fFitTotal->FixParameter(11,0.1);
387 }
388 else
389 {
390 fFitTotal->SetParameter(10,rpsiprime->Parameter(1)); fFitTotal->SetParLimits(10,3.5,3.8); // mean'
391 fFitTotal->SetParameter(11,rpsiprime->Parameter(2)); fFitTotal->SetParLimits(11,0.05,0.7); // sigma'
392 }
393
394 TFitResultPtr rpol1 = h.Fit("pol1",fitOption,"",1.5,2.5);
395 fFitTotal->SetParameter( 7,rpol1->Parameter(0));
396 fFitTotal->SetParameter( 8,rpol1->Parameter(1));
397
398 TFitResultPtr rexpo = h.Fit("expo",fitOption,"",4.5,7.0);
399 fFitTotal->SetParameter(12,rexpo->Parameter(0));
400 fFitTotal->SetParameter(13,rexpo->Parameter(1));
401
402
403 TFitResultPtr r = h.Fit(fFitTotal,fitOption,"",1.5,7);
404
405 TF1* signal = new TF1("signal","gaus",2,6);
406 signal->SetParameters(fFitTotal->GetParameter(0),
407 fFitTotal->GetParameter(3),
408 fFitTotal->GetParameter(4));
409
410 TF1* signalPrime = new TF1("signalPrime","gaus",2,6);
411 signalPrime->SetParameters(fFitTotal->GetParameter(9),
412 fFitTotal->GetParameter(10),
413 fFitTotal->GetParameter(11));
414
415 Double_t gausParameters[3];
416 Double_t covarianceMatrix[3][3];
417 Double_t gausParametersPrime[3];
418 Double_t covarianceMatrixPrime[3][3];
419
420 covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
421 covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
422 covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
423 covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
424 covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
425
426 for ( int iy = 1; iy < 3; ++iy )
427 {
428 for ( int ix = 1; ix < 3; ++ix )
429 {
430 covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
431 }
432 }
433
434 gausParameters[0] = fFitTotal->GetParameter(0);
435 gausParameters[1] = fFitTotal->GetParameter(3);
436 gausParameters[2] = fFitTotal->GetParameter(4);
437
438 gausParametersPrime[0] = fFitTotal->GetParameter(9);
439 gausParametersPrime[1] = fFitTotal->GetParameter(10);
440 gausParametersPrime[2] = fFitTotal->GetParameter(11);
441
442 covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(9,9);
443 covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(10,9);
444 covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(11,9);
445 covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(9,10);
446 covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(9,11);
447
448 for ( int iy = 1; iy < 3; ++iy )
449 {
450 for ( int ix = 1; ix < 3; ++ix )
451 {
452 covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
453 }
454 }
455
456 double n = signal->Integral(2,6)/h.GetBinWidth(10);
457 double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
458 Set("NofJpsi",n,nerr);
459 Set("MeanJpsi",fFitTotal->GetParameter(3),fFitTotal->GetParError(3));
460 Set("SigmaJpsi",fFitTotal->GetParameter(4),fFitTotal->GetParError(4));
461
462 double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
463 double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
464 Set("NofJpsiPrime",nprime,nerrprime);
465 Set("MeanJpsiPrime",fFitTotal->GetParameter(10),fFitTotal->GetParError(10));
466 Set("SigmaJpsiPrime",fFitTotal->GetParameter(11),fFitTotal->GetParError(11));
467}
468
469//_____________________________________________________________________________
470void AliAnalysisMuMu::Result::FitJpsiJpsiPrimeCB(TH1& h)
471{
472 std::cout << "Fit with jpsi + jpsiprime' (CB) " << std::endl;
473
474 const Double_t xmin(1.5);
475 const Double_t xmax(8.0);
476
477 fFitTotal = new TF1("fFitTotal",funcJpsiJpsiPrime,xmin,xmax,14);
478
479// Double_t N = par[0];
480// Double_t alpha = par[1];
481// Double_t n = par[2];
482// Double_t mean = par[3];
483// Double_t sigma = par[4];
484
485 fFitTotal->SetParameter( 0,1); // N
486 fFitTotal->FixParameter( 1,0.936); // alpha
487 fFitTotal->FixParameter( 2,4.44); // n
488 fFitTotal->SetParameter( 3,3.1); fFitTotal->SetParLimits(3,3.0,3.2); // mean
489 fFitTotal->SetParameter( 4,0.07); fFitTotal->SetParLimits(4,0.02,1); // sigma
490
491 fFitTotal->SetParameter( 5,0.01); // N'
492 fFitTotal->FixParameter( 6,0.936); // alpha'
493 fFitTotal->FixParameter( 7,4.44); // n'
494 fFitTotal->SetParameter( 8,3.7); fFitTotal->SetParLimits(8,3.5,3.8); // mean'
495 fFitTotal->SetParameter( 9,0.1); fFitTotal->SetParLimits(9,0.02,1.0); // sigma'
496
497 fFitTotal->SetParameter(10,h.GetMaximum());
498 fFitTotal->SetParameter(11,-1);
499
500 fFitTotal->SetParameter(12,h.GetMaximum()/100);
501 fFitTotal->SetParameter(13,-1);
502
503 TFitResultPtr r = h.Fit(fFitTotal,"SQBI","",1.5,6);
504
505// for ( int ix = 0; ix < fFitTotal->GetNpar(); ++ix )
506// {
507// for ( int iy = 0; iy < fFitTotal->GetNpar(); ++iy )
508// {
509// std::cout << Form("COV(%d,%d)=%e ",ix,iy,r->GetCovarianceMatrix()(ix,iy));
510// }
511// std::cout << std::endl;
512// }
513
514
515 TF1* signal = new TF1("signal","gaus",2,8);
516
517 signal->SetParameters(fFitTotal->GetParameter(0),
518 fFitTotal->GetParameter(3),
519 fFitTotal->GetParameter(4));
520
521 TF1* signalPrime = new TF1("signalPrime","gaus",2,8);
522
523 signalPrime->SetParameters(fFitTotal->GetParameter(0),
524 fFitTotal->GetParameter(8),
525 fFitTotal->GetParameter(9));
526
527 Double_t gausParameters[3];
528 Double_t gausParametersPrime[3];
529 Double_t covarianceMatrix[3][3];
530 Double_t covarianceMatrixPrime[3][3];
531
532 gausParameters[0] = fFitTotal->GetParameter(0);
533 gausParameters[1] = fFitTotal->GetParameter(3);
534 gausParameters[2] = fFitTotal->GetParameter(4);
535
536 covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
537 covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
538 covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
539 covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
540 covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
541
542 for ( int iy = 1; iy < 3; ++iy )
543 {
544 for ( int ix = 1; ix < 3; ++ix )
545 {
546 covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
547 }
548 }
549
550 gausParametersPrime[0] = fFitTotal->GetParameter(0);
551 gausParametersPrime[1] = fFitTotal->GetParameter(8);
552 gausParametersPrime[2] = fFitTotal->GetParameter(9);
553
554 covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(0,0);
555 covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(8,0);
556 covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(9,0);
557 covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(0,8);
558 covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(0,9);
559
560 for ( int iy = 1; iy < 3; ++iy )
561 {
562 for ( int ix = 1; ix < 3; ++ix )
563 {
564 covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
565 }
566 }
567
568 double n = signal->Integral(2,6)/h.GetBinWidth(10);
569 double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
570
571 Set("NofJpsi",n,nerr);
572 Set("MeanJpsi",fFitTotal->GetParameter(3),fFitTotal->GetParError(3));
573 Set("SigmaJpsi",fFitTotal->GetParameter(4),fFitTotal->GetParError(4));
574
575 double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
576 double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
577
578 Set("NofJpsiPrime",nprime,nerrprime);
579 Set("MeanJpsiPrime",fFitTotal->GetParameter(8),fFitTotal->GetParError(8));
580 Set("SigmaJpsiPrime",fFitTotal->GetParameter(9),fFitTotal->GetParError(9));
581
582}
583
584//_____________________________________________________________________________
585void AliAnalysisMuMu::Result::FitJpsiGCBE(TH1& h)
586{
587 std::cout << "Fit with jpsi alone (gaus + CB + expo)" << std::endl;
588
589 const Double_t xmin(1.0);
590 const Double_t xmax(8.0);
591
592 fFitTotal = new TF1("fFitTotal",funcJpsiGCBE,xmin,xmax,10);
593 fFitTotal->SetParNames("cste","x0","sigma0","N","alpha","n","mean","sigma","expocste","exposlope");
594
595 fFitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
596
597 const Double_t cbalpha(0.98);
598 const Double_t cbn(5.2);
599
600 fFitTotal->FixParameter(4,cbalpha);
601 fFitTotal->FixParameter(5,cbn);
602
603 fFitTotal->SetParLimits(6,2.8,3.2); // mean
604 fFitTotal->SetParLimits(7,0.02,0.3); // sigma
605
606 TF1* fg = new TF1("fg","gaus",xmin,xmax);
607
608 h.Fit(fg,"","",0.75,3.0);
609
610 fFitTotal->SetParameter(0,fg->GetParameter(0));
611 fFitTotal->SetParameter(1,fg->GetParameter(1));
612 fFitTotal->SetParameter(2,fg->GetParameter(2));
613
614 TF1* fexpo = new TF1("expo","expo",xmin,xmax);
615
616 h.Fit(fexpo,"","",3.5,5);
617
618 fFitTotal->SetParameter(8,fexpo->GetParameter(0));
619 fFitTotal->SetParameter(9,fexpo->GetParameter(1));
620
621 fFitTotal->SetParameter(3,h.GetMaximum()),
622 fFitTotal->SetParameter(4,cbalpha);
623 fFitTotal->SetParameter(5,cbn);
624 fFitTotal->SetParameter(6,3.15);
625 fFitTotal->SetParameter(7,0.1);
626
627 const char* fitOption = "SI+";
628
629 TFitResultPtr r = h.Fit(fFitTotal,fitOption,"",2,5);
630
631 Set("MeanJpsi",fFitTotal->GetParameter(6),fFitTotal->GetParError(6));
632 Set("SigmaJpsi",fFitTotal->GetParameter(7),fFitTotal->GetParError(7));
633
634 double m = MeanJpsi();
635 double s = SigmaJpsi();
636 double n = 3.0;
637
638 TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
639 fcb->SetParameters(fFitTotal->GetParameter(3),
640 fFitTotal->GetParameter(4),
641 fFitTotal->GetParameter(5),
642 fFitTotal->GetParameter(6),
643 fFitTotal->GetParameter(7));
644
645 fcb->SetLineColor(6);
646 fcb->SetNpx(100);
647 TLine* l1 = new TLine(m-n*s,0,m-n*s,fFitTotal->GetParameter(3));
648 TLine* l2 = new TLine(m+n*s,0,m+n*s,fFitTotal->GetParameter(3));
649 l1->SetLineColor(6);
650 l2->SetLineColor(6);
651 h.GetListOfFunctions()->Add(fcb);
652 h.GetListOfFunctions()->Add(l1);
653 h.GetListOfFunctions()->Add(l2);
654
655
656 Double_t cbParameters[5];
657 Double_t covarianceMatrix[5][5];
658
659 cbParameters[0] = fFitTotal->GetParameter(3);
660 cbParameters[1] = fFitTotal->GetParameter(4);
661 cbParameters[2] = fFitTotal->GetParameter(5);
662 cbParameters[3] = fFitTotal->GetParameter(6);
663 cbParameters[4] = fFitTotal->GetParameter(7);
664
665 for ( int iy = 0; iy < 5; ++iy )
666 {
667 for ( int ix = 0; ix < 5; ++ix )
668 {
669 covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+3,iy+3);
670 }
671 }
672
673 double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
674
675 double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
676
677 Set("NofJpsi",njpsi,nerr);
678}
679
680//_____________________________________________________________________________
681void AliAnalysisMuMu::Result::FitJpsiPCBE(TH1& h)
682{
683 std::cout << "Fit with jpsi alone (pol2 + CB + expo)" << std::endl;
684
685 const Double_t xmin(2.0);
686 const Double_t xmax(5.0);
687
688 fFitTotal = new TF1("fFitTotal",funcJpsiPCBE,xmin,xmax,10);
689 fFitTotal->SetParNames("p0","p1","p2","N","alpha","n","mean","sigma","expocste","exposlope");
690
691 fFitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
692
693 const Double_t cbalpha(0.98);
694 const Double_t cbn(5.2);
695
696 fFitTotal->FixParameter(4,cbalpha);
697 fFitTotal->FixParameter(5,cbn);
698
699 fFitTotal->SetParLimits(6,2,4); // mean
700 fFitTotal->SetParLimits(7,0.05,0.2); // sigma
701
702 TF1* fpol2 = new TF1("pol2","pol2",xmin,xmax);
703
704 h.Fit(fpol2,"+","",2,2.8);
705
706 fFitTotal->SetParameter(0,fpol2->GetParameter(0));
707 fFitTotal->SetParameter(1,fpol2->GetParameter(1));
708 fFitTotal->SetParameter(2,fpol2->GetParameter(2));
709
710 TF1* fexpo = new TF1("expo","expo",xmin,xmax);
711
712 h.Fit(fexpo,"+","",3.5,4.5);
713
714 fFitTotal->SetParameter(8,fexpo->GetParameter(0));
715 fFitTotal->SetParameter(9,fexpo->GetParameter(1));
716
717 fFitTotal->SetParameter(3,h.GetMaximum()),
718 fFitTotal->SetParameter(4,cbalpha);
719 fFitTotal->SetParameter(5,cbn);
720 fFitTotal->SetParameter(6,3.15);
721 fFitTotal->SetParameter(7,0.1);
722
723 h.Fit(fFitTotal,"+","",2.5,5);
724
725 Set("MeanJpsi",fFitTotal->GetParameter(6),fFitTotal->GetParError(6));
726 Set("SigmaJpsi",fFitTotal->GetParameter(7),fFitTotal->GetParError(7));
727
728 double m = MeanJpsi();
729 double s = SigmaJpsi();
730 double n = 2.0;
731
732 TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
733 fcb->SetParameters(fFitTotal->GetParameter(3),
734 fFitTotal->GetParameter(4),
735 fFitTotal->GetParameter(5),
736 fFitTotal->GetParameter(6),
737 fFitTotal->GetParameter(7));
738
739 fcb->SetLineColor(6);
740 fcb->SetNpx(100);
741 TLine* l1 = new TLine(m-n*s,0,m-n*s,fFitTotal->GetParameter(3));
742 TLine* l2 = new TLine(m+n*s,0,m+n*s,fFitTotal->GetParameter(3));
743 l1->SetLineColor(6);
744 l2->SetLineColor(6);
745 h.GetListOfFunctions()->Add(fcb);
746 h.GetListOfFunctions()->Add(l1);
747 h.GetListOfFunctions()->Add(l2);
748
749
750 Set("NofJpsi",fFitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fFitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
751
752 // Set("NofJpsi",fFitTotal->Integral(0,10)/h.GetBinWidth(1),fFitTotal->IntegralError(0,10)/h.GetBinWidth(1));
753
754}
755
756//_____________________________________________________________________________
757void AliAnalysisMuMu::Result::FitJpsiCBE(TH1& h)
758{
759 std::cout << "Fit with jpsi alone" << std::endl;
760
761 const Double_t xmin(1.5);
762 const Double_t xmax(8.0);
763
764 fFitTotal = new TF1("fFitTotal",funcJpsiCBE,xmin,xmax,7);
765 fFitTotal->SetParNames("N","alpha","n","mean","sigma","expocste","exposlope");
766
767// fFitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3,1,0);
768
769 fFitTotal->SetParameters(1,1.15,3.6,3.0,0.07,1,-1);
770
771 fFitTotal->SetParLimits(0,0,h.GetMaximum()); // N
772// fFitTotal->SetParLimits(1,0.1,2); // alpha
773 fFitTotal->FixParameter(1,0.98);
774// fFitTotal->SetParLimits(2,0.01,5); // n
775 fFitTotal->FixParameter(2,5.2);
776 fFitTotal->SetParLimits(3,2.8,3.5); // mean
777 fFitTotal->SetParLimits(4,0.05,0.2); // sigma
778
779 TF1* fexpo = new TF1("expo","expo",xmin,xmax);
780
781 h.Fit(fexpo,"+","",2,3);
782
783 fFitTotal->SetParameter(5,fexpo->GetParameter(0));
784 fFitTotal->SetParameter(6,fexpo->GetParameter(1));
785
786 h.Fit(fFitTotal,"+","",2,5);
787
788
789 Set("MeanJpsi",fFitTotal->GetParameter(3),fFitTotal->GetParError(3));
790 Set("SigmaJpsi",fFitTotal->GetParameter(4),fFitTotal->GetParError(4));
791
792 double m = MeanJpsi();
793 double s = SigmaJpsi();
794 double n = 2.0;
795
796 TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
797 fcb->SetParameters(fFitTotal->GetParameter(0),
798 fFitTotal->GetParameter(1),
799 fFitTotal->GetParameter(2),
800 fFitTotal->GetParameter(3),
801 fFitTotal->GetParameter(4));
802
803 fcb->SetLineColor(6);
804 fcb->SetNpx(1000);
805 TLine* l1 = new TLine(m-n*s,0,m-n*s,fFitTotal->GetParameter(0));
806 TLine* l2 = new TLine(m+n*s,0,m+n*s,fFitTotal->GetParameter(0));
807 l1->SetLineColor(6);
808 l2->SetLineColor(6);
809 h.GetListOfFunctions()->Add(fcb);
810 h.GetListOfFunctions()->Add(l1);
811 h.GetListOfFunctions()->Add(l2);
812
813
814 Set("NofJpsi",fFitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fFitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
815
816 // Set("NofJpsi",fFitTotal->Integral(0,10)/h.GetBinWidth(1),fFitTotal->IntegralError(0,10)/h.GetBinWidth(1));
817
818}
819
820//_____________________________________________________________________________
821void AliAnalysisMuMu::Result::FitJpsiECBE(TH1& h)
822{
823 std::cout << "Fit with jpsi alone (expo + CB + expo)" << std::endl;
824
825 const Double_t xmin(1.5);
826 const Double_t xmax(8.0);
827
828 fFitTotal = new TF1("fFitTotal",funcJpsiECBE,xmin,xmax,9);
829 fFitTotal->SetParNames("e0","s0","N","alpha","n","mean","sigma","e1","s1");
830
831 fFitTotal->SetParameters(1,-1,1,1.15,3.6,3.2,0.06,-1);
832
833 fFitTotal->SetParLimits(0,0,h.GetMaximum()*2);
834
835 fFitTotal->FixParameter(3,0.98); // alpha
836 fFitTotal->FixParameter(4,5.2); // n
837 fFitTotal->SetParLimits(5,2.8,3.5); // mean
838 fFitTotal->SetParLimits(6,0.05,0.2); // sigma
839
840 TF1* fexpo1 = new TF1("expo1","expo",xmin,xmax);
841 TF1* fexpo2 = new TF1("expo2","expo",xmin,xmax);
842
843 h.Fit(fexpo1,"","",1.5,3);
844
845 fFitTotal->SetParameter(0,fexpo1->GetParameter(0));
846 fFitTotal->SetParameter(1,fexpo1->GetParameter(1));
847
848 h.Fit(fexpo2,"","",3.5,5.0);
849
850 fFitTotal->SetParameter(7,fexpo2->GetParameter(0));
851 fFitTotal->SetParameter(8,fexpo2->GetParameter(1));
852
853 const char* fitOption = "SI+";
854
855 TFitResultPtr r = h.Fit(fFitTotal,fitOption,"",2,5);
856
857 Set("MeanJpsi",fFitTotal->GetParameter(5),fFitTotal->GetParError(5));
858 Set("SigmaJpsi",fFitTotal->GetParameter(6),fFitTotal->GetParError(6));
859
860 double m = MeanJpsi();
861 double s = SigmaJpsi();
862 double n = 3.0;
863
864 TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
865 fcb->SetParameters(fFitTotal->GetParameter(2),
866 fFitTotal->GetParameter(3),
867 fFitTotal->GetParameter(4),
868 fFitTotal->GetParameter(5),
869 fFitTotal->GetParameter(6));
870
871 fcb->SetParError(0,fFitTotal->GetParError(2));
872 fcb->SetParError(1,fFitTotal->GetParError(3));
873 fcb->SetParError(2,fFitTotal->GetParError(4));
874 fcb->SetParError(3,fFitTotal->GetParError(5));
875 fcb->SetParError(4,fFitTotal->GetParError(6));
876
877 fcb->SetLineColor(6);
878 fcb->SetNpx(1000);
879 TLine* l1 = new TLine(m-n*s,0,m-n*s,fFitTotal->GetParameter(2));
880 TLine* l2 = new TLine(m+n*s,0,m+n*s,fFitTotal->GetParameter(2));
881 l1->SetLineColor(6);
882 l2->SetLineColor(6);
883 h.GetListOfFunctions()->Add(fcb);
884 h.GetListOfFunctions()->Add(l1);
885 h.GetListOfFunctions()->Add(l2);
886
887
888 Double_t cbParameters[5];
889 Double_t covarianceMatrix[5][5];
890
891 cbParameters[0] = fFitTotal->GetParameter(2);
892 cbParameters[1] = fFitTotal->GetParameter(3);
893 cbParameters[2] = fFitTotal->GetParameter(4);
894 cbParameters[3] = fFitTotal->GetParameter(5);
895 cbParameters[4] = fFitTotal->GetParameter(6);
896
897 for ( int iy = 0; iy < 5; ++iy )
898 {
899 for ( int ix = 0; ix < 5; ++ix )
900 {
901 covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
902 }
903 }
904
905
906 double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
907
908 double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
909
910
911 Set("NofJpsi",njpsi,nerr);
912}
913
914//_____________________________________________________________________________
915void AliAnalysisMuMu::Result::FitJpsi(TH1& h)
916{
917 std::cout << "Fit with jpsi alone" << std::endl;
918
919 const Double_t xmin(1.5);
920 const Double_t xmax(8.0);
921
922 fFitTotal = new TF1("fFitTotal",funcCB2,xmin,xmax,7);
923 fFitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
924 fFitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3);
925 fFitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
926 fFitTotal->SetParLimits(1,0,10); // alpha
927 fFitTotal->SetParLimits(2,1,10); // n
928 fFitTotal->SetParLimits(3,1,4); // mean
929 fFitTotal->SetParLimits(4,0.01,1); // sigma
930 fFitTotal->SetParLimits(5,0,10); // alpha
931 fFitTotal->SetParLimits(6,1,10); // n
932
933 h.Fit(fFitTotal,"+","",2,5);
934
935
936 Set("MeanJpsi",fFitTotal->GetParameter(3),fFitTotal->GetParError(3));
937 Set("SigmaJpsi",fFitTotal->GetParameter(4),fFitTotal->GetParError(4));
938
939 double m = MeanJpsi();
940 double s = SigmaJpsi();
941 double n = 2.0;
942
943 Set("NofJpsi",fFitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fFitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
944
945// Set("NofJpsi",fFitTotal->Integral(0,10)/h.GetBinWidth(1),fFitTotal->IntegralError(0,10)/h.GetBinWidth(1));
946
947}
948
949//_____________________________________________________________________________
950void AliAnalysisMuMu::Result::FitUpsilon(TH1& h)
951{
952 std::cout << "Fit with upsilon alone" << std::endl;
953
954 const Double_t xmin(6.0);
955 const Double_t xmax(12.0);
956
957 fFitTotal = new TF1("fFitTotal",funcCB2,xmin,xmax,7);
958 fFitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
959 fFitTotal->SetParameters(h.GetMaximum(),1,5,9.46,0.2,1.5,3);
960 fFitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
961 fFitTotal->SetParLimits(1,0,10); // alpha
962 fFitTotal->SetParLimits(2,1,10); // n
963 fFitTotal->SetParLimits(3,8,12); // mean
964 fFitTotal->SetParLimits(4,0.01,1); // sigma
965 fFitTotal->SetParLimits(5,0,10); // alpha
966 fFitTotal->SetParLimits(6,1,10); // n
967
968 h.Fit(fFitTotal,"+","",6,12);
969
970
971 Set("MeanUpsilon",fFitTotal->GetParameter(3),fFitTotal->GetParError(3));
972 Set("SigmaUpsilon",fFitTotal->GetParameter(4),fFitTotal->GetParError(4));
973
974 double m = MeanUpsilon();
975 double s = SigmaUpsilon();
976 double n = 2.0;
977
978 Set("NofUpsilon",fFitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fFitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
979}
980
981//_____________________________________________________________________________
982void AliAnalysisMuMu::Result::Fit(Int_t nrebin)
983{
984 // new TCanvas;
985
986 static int n(0);
987
988 fRebin = nrebin;
989
990// TH2* h2 = static_cast<TH2*>(fHC->Histo("MinvUSPt"));
991// int binptmin = h2->GetXaxis()->FindBin(2.0);
992// TH1* hminv = static_cast<TH1*>(h2->ProjectionY("minv0",binptmin,h2->GetXaxis()->GetNbins()));
993
994 TH1* hminv = fHC->Histo("MinvUSPt:py");
995
2754fd07 996 if (!hminv || hminv->GetEntries()<100 ) return;
2f331ac9 997
998 fMinv = static_cast<TH1*>(hminv->Clone(Form("minv%d",n)));
999 fMinv->Rebin(fRebin);
1000 fMinv->SetDirectory(0);
1001// fMinv->SetStats(false);
1002 fMinv->SetTitle(TriggerClass());
1003
1004 if ( ( fFitType & kJpsi ) &&
1005 ( fFitType & kJpsiPrime ) )
1006 {
1007 FitJpsiJpsiPrimeCustom(*fMinv);
1008// FitJpsiJpsiPrimeCB(*fMinv);
1009 }
1010 else if ( ( fFitType & kJpsi ) && ( fFitType & kPbPb2011 ) )
1011 {
1012 if ( fFitType & kMatchAny )
1013 {
1014 FitJpsiECBE(*fMinv); // for MATCH
1015 }
1016 else
1017 {
1018 // FitJpsiPCBE(*fMinv);
1019 FitJpsiGCBE(*fMinv); // for MATCHLOW
1020 }
1021 fMinv->GetXaxis()->SetRangeUser(2.0,5.0);
1022 fMinv->SetMarkerStyle(kFullDotMedium);
1023 fMinv->SetMarkerColor(kBlue);
1024 fMinv->SetLineColor(kBlue);
1025 fMinv->SetDrawOption("EP");
1026 }
1027 else
1028 {
1029 int bminjpsi = fMinv->GetXaxis()->FindBin(2);
1030 int bmaxjpsi = fMinv->GetXaxis()->FindBin(4);
1031
1032 int bminupsilon = fMinv->GetXaxis()->FindBin(8);
1033 int bmaxupsilon = fMinv->GetXaxis()->FindBin(12);
1034
1035 if ( fMinv->Integral(bminjpsi,bmaxjpsi) > 100.0 ) FitJpsi(*fMinv);
1036 if ( fMinv->Integral(bminupsilon,bmaxupsilon) > 100.0 ) FitUpsilon(*fMinv);
1037 }
1038
1039 TH2* hpp = static_cast<TH2*>(fHC->Histo("MinvPPPt"));
1040 TH2* hmm = static_cast<TH2*>(fHC->Histo("MinvMMPt"));
1041
1042 if ( hpp && hmm )
1043 {
1044 TH2* htmp = static_cast<TH2*>(hpp->Clone("htmp"));
1045 htmp->SetDirectory(0);
1046 htmp->Add(hmm);
1047 htmp->Scale(0.5);
1048
1049 fMinvLS = htmp->ProjectionY();
1050 fMinvLS->SetDirectory(0);
1051 fMinvLS->Rebin(fRebin);
1052
1053 delete htmp;
1054 }
1055 else
1056 {
1057 AliDebug(1,Form("hpp=%p hmm=%p",hpp,hmm)); // this might happen for small runs and
1058 // strict cuts on pairs
1059 }
1060
2754fd07 1061// Print();
2f331ac9 1062}
1063
1064//_____________________________________________________________________________
1065TMap* AliAnalysisMuMu::Result::Map() const
1066{
1067 if (!fMap)
1068 {
1069 fMap = new TMap;
1070 fMap->SetOwnerKeyValue(kTRUE,kTRUE);
1071 }
1072 return fMap;
1073}
1074
1075//_____________________________________________________________________________
1076void AliAnalysisMuMu::Result::Set(const char* name, double value, double error)
1077{
1078 TObjArray* p = static_cast<TObjArray*>(Map()->GetValue(name));
1079 if (!p)
1080 {
1081 TParameter<double>* v = new TParameter<double>(name,value);
1082 TParameter<double>* e = new TParameter<double>(name,error);
1083 p = new TObjArray;
1084 p->SetOwner(kTRUE);
1085 p->Add(v);
1086 p->Add(e);
1087 Map()->Add(new TObjString(name),p);
1088 }
1089 else
1090 {
1091 TParameter<double>* v = static_cast<TParameter<double>*>(p->At(0));
1092 TParameter<double>* e = static_cast<TParameter<double>*>(p->At(1));
1093 v->SetVal(value);
1094 e->SetVal(error);
1095 }
1096}
1097
1098//_____________________________________________________________________________
1099Bool_t AliAnalysisMuMu::Result::HasValue(const char* name) const
1100{
1101 return ( Map()->GetValue(name) != 0x0 );
1102}
1103
1104//_____________________________________________________________________________
1105Double_t AliAnalysisMuMu::Result::GetValue(const char* name) const
1106{
1107 TObjArray* p = static_cast<TObjArray*>(Map()->GetValue(name));
1108 if (p)
1109 {
1110 TParameter<double>* val = static_cast<TParameter<double>*>(p->At(0));
1111 return val->GetVal();
1112 }
1113 return 0.0;
1114}
1115
1116//_____________________________________________________________________________
1117Double_t AliAnalysisMuMu::Result::GetError(const char* name) const
1118{
1119 TObjArray* p = static_cast<TObjArray*>(Map()->GetValue(name));
1120 if (p)
1121 {
1122 TParameter<double>* val = static_cast<TParameter<double>*>(p->At(1));
1123 return val->GetVal();
1124 }
1125 return 0.0;
1126}
1127
1128//_____________________________________________________________________________
1129void AliAnalysisMuMu::Result::Print(Option_t* /*opt*/) const
1130{
2754fd07 1131 std::cout << Form("%20s - %s %s %s - NRUNS %d - NTRIGGER %10d",
2f331ac9 1132 GetName(),
2754fd07 1133 EventSelection(),
1134 PairSelection(),
1135 CentralitySelection(),
1136 NofRuns(),
1137 NofTriggers()) << std::endl;
2f331ac9 1138
1139 if (NofJpsi())
1140 {
1141 std::cout << Form("\t\tNjpsi %7.2f +- %5.2f \n\t\tRatio (x10^4) %7.2f +- %5.2f CHI2/NDF %5.2f"
1142 "\n\t\tPEAK %7.2f +- %5.2f Gev/c^2 SIGMA %7.2f +- %5.2f MeV/c^2",
1143 NofJpsi(),ErrorOnNofJpsi(),
1144 NofJpsi()*1E4/fNofTriggers,
1145 ErrorOnNofJpsi()*1E4/fNofTriggers,
1146 (fFitTotal ? fFitTotal->GetChisquare()/fFitTotal->GetNDF() : 0),
1147 MeanJpsi(),ErrorOnMeanJpsi(),
1148 1E3*SigmaJpsi(),1E3*ErrorOnSigmaJpsi()) << std::endl;
1149
1150 }
1151
1152 if (NofJpsiPrime())
1153 {
1154 std::cout << Form("\t\tNjpsiPrime %7.2f +- %5.2f \n\t\tRatio (x10^4) %7.2f +- %5.2f CHI2/NDF %5.2f"
1155 "\n\t\tPEAK %7.2f +- %5.2f Gev/c^2 SIGMA %7.2f +- %5.2f MeV/c^2",
1156 NofJpsiPrime(),ErrorOnNofJpsiPrime(),
1157 NofJpsiPrime()*1E4/fNofTriggers,
1158 ErrorOnNofJpsiPrime()*1E4/fNofTriggers,
1159 (fFitTotal ? fFitTotal->GetChisquare()/fFitTotal->GetNDF() : 0),
1160 MeanJpsiPrime(),ErrorOnMeanJpsiPrime(),
1161 1E3*SigmaJpsiPrime(),1E3*ErrorOnSigmaJpsiPrime()) << std::endl;
1162
1163 }
1164
1165 if ( NofUpsilon() )
1166 {
1167 std::cout << Form("\t\tNupsilon %7.2f +- %5.2f \n\t\tRatio (x10^4) %7.2f +- %5.2f CHI2/NDF %5.2f"
1168 "\n\t\tPEAK %7.2f +- %5.2f Gev/c^2 SIGMA %7.2f +- %5.2f MeV/c^2",
1169 NofUpsilon(),ErrorOnNofUpsilon(),
1170 NofUpsilon()*1E4/fNofTriggers,
1171 ErrorOnNofUpsilon()*1E4/fNofTriggers,
1172 (fFitTotal ? fFitTotal->GetChisquare()/fFitTotal->GetNDF() : 0),
1173 MeanUpsilon(),ErrorOnMeanUpsilon(),
1174 1E3*SigmaUpsilon(),1E3*ErrorOnSigmaUpsilon()) << std::endl;
1175 }
1176}
1177
1178
1179//_____________________________________________________________________________
2754fd07 1180TString FindTrigger(const AliHistogramCollection& hc,
2f331ac9 1181 const char* base,
1182 const char* selection,
1183 const char* paircut,
1184 const char* centrality)
1185{
1186 /// find the trigger containing the MinvPt histograms
1187
1188 std::vector<std::string> trigger2test;
1189
1190 // trigger2test.push_back(Form("%s5-B-NOPF-ALLNOTRD",base));
1191 // trigger2test.push_back(Form("%s1-B-NOPF-ALLNOTRD",base));
1192 // trigger2test.push_back(Form("%s1B-ABCE-NOPF-MUON",base));
1193 if ( TString(base).Contains("||") || TString(base).Contains("-") )
1194 {
1195 trigger2test.push_back(base);
1196 }
1197 else
1198 {
1199 trigger2test.push_back(Form("%s-B-NOPF-ALLNOTRD",base));
1200 trigger2test.push_back(Form("%s-B-NOPF-MUON",base));
1201 trigger2test.push_back(Form("%s-S-NOPF-ALLNOTRD",base));
1202 trigger2test.push_back(Form("%s-S-NOPF-MUON",base));
1203 }
1204 trigger2test.push_back("ANY");
1205
1206 for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
1207 {
1208 std::string trigger = trigger2test[i];
1209
1210 if ( hc.Histo(selection,trigger.c_str(),centrality,paircut,"MinvUSPt") )
1211 {
1212 return trigger.c_str();
1213 }
1214 }
1215
1216// AliWarningGeneral("FindTrigger",Form("DID NOT FIND TRIGGER base=%s selection=%s paircut=%s centrality=%s",
1217// base,selection,paircut,centrality));
1218// for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
1219// {
1220// AliWarningGeneral("FindTrigger",Form("tested trigger = %s",trigger2test[i].c_str()));
1221// }
1222 return "";
1223}
1224
1225//_____________________________________________________________________________
1226//_____________________________________________________________________________
1227//_____________________________________________________________________________
1228//_____________________________________________________________________________
1229//_____________________________________________________________________________
1230
1231TString AliAnalysisMuMu::fgOCDBPath("raw://");
1232
2754fd07 1233TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-MUON,CMUL8-S-NOPF-MUON,CMUL7-B-NOPF-ALLNOTRD,CMUL7-B-NOPF-MUON,CMUU7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-MUON,CPBI1MUL-B-NOPF-MUON");
1234
1235TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON,CMSL7-S-NOPF-ALLNOTRD,CMSL8-S-NOPF-MUON,CMSL8-S-NOPF-ALLNOTRD,CMSL7-B-NOPF-MUON,CMUS1-B-NOPF-MUON,CMUS7-B-NOPF-MUON");
1236
1237TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
1238
1239TString AliAnalysisMuMu::fgDefaultEventSelectionList("ALL");
1240
1241TString AliAnalysisMuMu::fgDefaultPairSelectionList("pMATCHLOWRABSDCABOTH");
1242
2f331ac9 1243//_____________________________________________________________________________
1244AliAnalysisMuMu::AliAnalysisMuMu(const char* filename) : TObject(),
1245fFilename(filename),
1246fHistogramCollection(0x0),
1247fCounterCollection(0x0),
2754fd07 1248fDimuonTriggers(fgDefaultDimuonTriggers),
1249fMuonTriggers(fgDefaultMuonTriggers),
1250fMinbiasTriggers(fgDefaultMinbiasTriggers),
1251fEventSelectionList(fgDefaultEventSelectionList),
1252fPairSelectionList(fgDefaultPairSelectionList)
2f331ac9 1253{
1254 // ctor
1255
1256 GetCollections(fFilename,fHistogramCollection,fCounterCollection);
1257}
1258
1259//_____________________________________________________________________________
1260AliAnalysisMuMu::~AliAnalysisMuMu()
1261{
1262 // dtor
1263 gROOT->CloseFiles();
1264 delete fCounterCollection;
1265 delete fHistogramCollection;
1266}
1267
1268//_____________________________________________________________________________
1269void AliAnalysisMuMu::CentralityCheck(const char* filelist)
1270{
1271 // Check if we get correctly filled centrality
1272
1273 TObjArray* files = ReadFileList(filelist);
1274
1275 if (!files || files->IsEmpty() ) return;
1276
1277 TIter next(files);
1278 TObjString* str;
1279
1280 while ( ( str = static_cast<TObjString*>(next()) ) )
1281 {
1282 AliHistogramCollection* hc(0x0);
1283 AliCounterCollection* cc(0x0);
1284
1285 if (!GetCollections(str->String().Data(),hc,cc)) continue;
1286
1287 int run = RunNumberFromFileName(str->String().Data());
1288
1289 TH1* h = hc->Histo("/ALL/CPBI1MUL-B-NOPF-MUON/Centrality");
1290
1291 float percent(0);
1292
1293 if (h)
1294 {
1295 percent = 100*h->Integral(1,1)/h->Integral();
1296 }
1297
1298 std::cout << Form("RUN %09d PERCENT %7.2f",run,percent) << std::endl;
1299
1300 delete hc;
1301 }
1302
1303 gROOT->CloseFiles();
1304
1305}
1306
1307//_____________________________________________________________________________
2754fd07 1308void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
1309 ULong64_t* totalNmb,
1310 ULong64_t* totalNmsl,
1311 ULong64_t* totalNmul)
2f331ac9 1312{
1313 // Report of some basic numbers, like number of MB and MUON triggers,
1314 // both before and after physics selection, and comparison with
1315 // the total number of such triggers (taken from the OCDB scalers)
1316 // if requested.
1317 //
1318 // filename is assumed to be a root filecontaining a list containing
1319 // an AliCounterCollection (or directly an AliCounterCollection)
1320 //
1321 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
1322 //
1323
1324 if (!fHistogramCollection || !fCounterCollection) return;
1325
1326 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
1327 TIter nextRun(runs);
1328
1329 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1330 TIter nextTrigger(triggers);
1331
1332 TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1333
2754fd07 1334 Bool_t doPS = (events->FindObject("PSALL") != 0x0);
1335
2f331ac9 1336 TObjString* srun;
1337 TObjString* strigger;
1338
2754fd07 1339 ULong64_t localNmb(0);
1340 ULong64_t localNmsl(0);
1341 ULong64_t localNmul(0);
1342
1343 if ( totalNmb) *totalNmb = 0;
1344 if ( totalNmsl) *totalNmsl = 0;
1345 if ( totalNmul ) *totalNmul = 0;
2f331ac9 1346
1347 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
1348 {
1349 std::cout << Form("RUN %09d ",srun->String().Atoi());
1350
1351 TString details;
1352 ULong64_t nmb(0);
1353 ULong64_t nmsl(0);
1354 ULong64_t nmul(0);
1355
1356 nextTrigger.Reset();
1357
2754fd07 1358 Int_t nofPS(0);
1359
2f331ac9 1360 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
1361 {
2754fd07 1362
1363 if ( !fgDefaultMinbiasTriggers.Contains(strigger->String().Data()) &&
1364 !fgDefaultMuonTriggers.Contains(strigger->String().Data()) &&
1365 !fgDefaultDimuonTriggers.Contains(strigger->String().Data()) ) continue;
1366
2f331ac9 1367 ULong64_t n = fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
1368 strigger->String().Data(),"ALL",srun->String().Atoi()));
1369
2754fd07 1370 details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
2f331ac9 1371
1372
2754fd07 1373 ULong64_t nps = fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
1374 strigger->String().Data(),"PSALL",srun->String().Atoi()));
1375
1376 if ( doPS )
2f331ac9 1377 {
2754fd07 1378 details += TString::Format(" PS %5.1f %%",nps*100.0/n);
2f331ac9 1379 }
1380
2754fd07 1381 if (nps)
1382 {
1383 ++nofPS;
1384 }
1385
2f331ac9 1386 if ( fMinbiasTriggers.Contains(strigger->String()) )
1387 {
1388 nmb += n;
2754fd07 1389 if ( totalNmb) (*totalNmb) += n;
1390 localNmb += n;
2f331ac9 1391 }
1392 else if ( fMuonTriggers.Contains(strigger->String()) )
1393 {
1394 nmsl += n;
2754fd07 1395 if ( totalNmsl) (*totalNmsl) += n;
1396 localNmsl += n;
2f331ac9 1397 }
1398 else if ( fDimuonTriggers.Contains(strigger->String()) )
1399 {
1400 nmul += n;
2754fd07 1401 if ( totalNmul ) (*totalNmul) += n;
1402 localNmul += n;
2f331ac9 1403 }
1404 }
1405
2754fd07 1406 std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
1407 nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
2f331ac9 1408
1409 if ( detailTriggers )
1410 {
1411 std::cout << details.Data();
1412 }
1413 std::cout << std::endl;
1414 }
1415
2754fd07 1416 if ( !totalNmul && !totalNmsl && !totalNmb )
1417 {
1418 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
1419 localNmb,localNmsl,localNmul) << std::endl;
1420 }
2f331ac9 1421
2754fd07 1422 delete runs;
2f331ac9 1423 delete triggers;
1424 delete events;
1425}
1426
2754fd07 1427//_____________________________________________________________________________
1428void AliAnalysisMuMu::BasicCountsEvolution(const char* filelist, Bool_t detailTriggers)
1429{
1430 // Report of some basic numbers, like number of MB and MUON triggers,
1431 // both before and after physics selection, and comparison with
1432 // the total number of such triggers (taken from the OCDB scalers)
1433 // if requested.
1434 //
1435 // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
1436 //
1437 // To change the list of (single muon, dimuon, MB) triggers, use
1438 // the SetDefault*TriggerList methods prior to call this one
1439 //
1440
1441 TObjArray* files = ReadFileList(filelist);
1442
1443 if (!files || files->IsEmpty() ) return;
1444
1445 TIter next(files);
1446 TObjString* str;
1447
1448 ULong64_t totalNmb(0);
1449 ULong64_t totalNmsl(0);
1450 ULong64_t totalNmul(0);
1451
1452 while ( ( str = static_cast<TObjString*>(next()) ) )
1453 {
1454 AliAnalysisMuMu m(str->String().Data());
1455
1456 ULong64_t nmb(0);
1457 ULong64_t nmsl(0);
1458 ULong64_t nmul(0);
1459
1460 m.BasicCounts(detailTriggers,&nmb,&nmsl,&nmul);
1461
1462 totalNmb += nmb;
1463 totalNmsl += nmsl;
1464 totalNmul += nmul;
1465 }
1466
1467 std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
1468 totalNmb,totalNmsl,totalNmul) << std::endl;
1469
1470}
1471
2f331ac9 1472//_____________________________________________________________________________
1473TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
1474 const char* backgroundresults)
1475{
1476 TFile* fjpsi = FileOpen(jpsiresults);
1477 TFile* fbck = FileOpen(backgroundresults);
1478
1479 if (!fjpsi || !fbck) return 0x0;
1480
1481 TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
1482
1483 std::vector<std::string> checks;
1484
1485 checks.push_back("muminus-CMUU7-B-NOPF-ALLNOTRD");
1486 checks.push_back("muplus-CMUU7-B-NOPF-ALLNOTRD");
1487 checks.push_back("muminus-CMUSH7-B-NOPF-MUON");
1488 checks.push_back("muplus-CMUSH7-B-NOPF-MUON");
1489
1490 if (!gjpsi) return 0x0;
1491
1492 TObjArray* a = new TObjArray;
1493 a->SetOwner(kTRUE);
1494
1495 for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
1496 {
1497
1498 TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
1499
1500 if (!gback) continue;
1501
1502 if ( gjpsi->GetN() != gback->GetN() )
1503 {
1504 AliErrorClass("graphs have different number of points !");
1505 continue;
1506 }
1507
1508 TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
1509
1510 for ( int i = 0; i < gjpsi->GetN(); ++i )
1511 {
1512 double r1,r2,y1,y2;
1513
1514 gjpsi->GetPoint(i,r1,y1);
1515 gback->GetPoint(i,r2,y2);
1516
1517 if ( r1 != r2 )
1518 {
1519 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
1520 continue;
1521 }
1522
1523 g->SetPoint(i,y2,y1);
1524 // g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
1525 }
1526
1527 g->SetMarkerStyle(25+j);
1528 g->SetMarkerSize(1.2);
1529 if (j==0)
1530 {
1531 g->Draw("ap");
1532 }
1533 else
1534 {
1535 g->Draw("p");
1536 }
1537 g->SetLineColor(j+1);
1538 g->SetMarkerColor(j+1);
1539 g->SetName(checks[j].c_str());
1540 a->AddLast(g);
1541 }
1542
1543 return a;
1544}
1545
1546//_____________________________________________________________________________
1547TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
1548 const char* simjpsiresults)
1549{
1550 TFile* freal = FileOpen(realjpsiresults);
1551 TFile* fsim = FileOpen(simjpsiresults);
1552
1553 if (!freal || !fsim) return 0x0;
1554
1555 TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
1556 TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
1557
1558 TObjArray* a = new TObjArray;
1559 a->SetOwner(kTRUE);
1560
1561 if ( greal->GetN() != gsim->GetN() )
1562 {
1563 AliErrorClass("graphs have different number of points !");
1564 return 0x0;
1565 }
1566
1567 TGraphErrors* g = new TGraphErrors(greal->GetN());
1568 TGraphErrors* gratio = new TGraphErrors(greal->GetN());
1569
1570 for ( int i = 0; i < greal->GetN(); ++i )
1571 {
1572 double r1,r2,y1,y2;
1573
1574 greal->GetPoint(i,r1,y1);
1575 gsim->GetPoint(i,r2,y2);
1576
1577 if ( r1 != r2 )
1578 {
1579 AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
1580 continue;
1581 }
1582
1583 double ratio(0.0);
1584
1585 if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
1586 {
1587 g->SetPoint(i,0,0);
1588 g->SetPointError(i,0,0);
1589 }
1590 else
1591 {
1592 g->SetPoint(i,y2,y1);
1593 g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
1594 ratio = y2/y1;
1595 }
1596 gratio->SetPoint(i,r1,ratio);
1597 }
1598
1599 g->SetMarkerStyle(25);
1600 g->SetMarkerSize(1.2);
1601
1602 new TCanvas;
1603
1604 g->Draw("ap");
1605
1606 g->SetLineColor(1);
1607 g->SetMarkerColor(1);
1608 g->SetName("jpsipercmuurealvssim");
1609
1610 new TCanvas;
1611
1612 greal->Draw("alp");
1613 gsim->SetLineColor(4);
1614
1615 gsim->Draw("lp");
1616
1617 new TCanvas;
1618 gratio->Draw("alp");
1619
1620 return g;
1621}
1622
1623//_____________________________________________________________________________
1624TObjArray* AliAnalysisMuMu::ComputeBackgroundEvolution(const char* filelist,
1625 const char* triggerList,
1626 const char* outputFile,
1627 const char* outputMode)
1628{
1629 // triggerList is a list of complete trigger names, separated by space
1630 // of the triggers to consider : only the first one found in the list
1631 // is used for each run.
1632
1633 TObjArray* files = ReadFileList(filelist);
1634
1635 if (!files || files->IsEmpty() ) return 0x0;
1636
1637 TIter next(files);
1638 TObjString* str;
1639
1640 const char* ps = "ALL";
1641 const char* centrality = "PP";
1642 const char* ts1 = "sMATCHLOWRABS";
1643 const char* ts2 = "sMATCHLOWRABSDCA";
1644
1645 std::map<std::string, std::vector<float> > runs;
1646 std::map<std::string, std::vector<float> > errruns;
1647 std::map<std::string, std::vector<float> > yplus,erryplus;
1648 std::map<std::string, std::vector<float> > yminus,erryminus;
1649
2754fd07 1650 TObjArray* triggers = TString(triggerList).Tokenize(",");
2f331ac9 1651
1652 TObjArray* a = new TObjArray;
1653 a->SetOwner(kTRUE);
1654
1655 Bool_t bothSigns(kFALSE);
1656
1657 while ( ( str = static_cast<TObjString*>(next()) ) )
1658 {
1659 AliInfoClass(str->String().Data());
1660
1661 AliHistogramCollection* hc(0x0);
1662 AliCounterCollection* cc(0x0);
1663
1664 if (!GetCollections(str->String().Data(),hc,cc)) continue;
1665
1666 TIter nextHisto(hc->CreateIterator());
1667 TH1* h;
1668 int nplus(0), nminus(0);
1669
1670 while ( ( h = static_cast<TH1*>(nextHisto()) ) )
1671 {
1672 if ( TString(h->GetName()).EndsWith("Plus") )
1673 {
1674 nplus++;
1675 }
1676 if ( TString(h->GetName()).EndsWith("Minus") )
1677 {
1678 nminus++;
1679 }
1680 }
1681
1682 if (nminus==nplus && nplus>0 )
1683 {
1684 bothSigns = kTRUE;
1685 }
1686
1687 AliInfoClass(Form("Both signs = %d",bothSigns));
1688
1689 TIter nextTrigger(triggers);
1690 TObjString* trigger;
1691 TH1* h1p(0x0);
1692 TH1* h1m(0x0);
1693 TH1* h2p(0x0);
1694 TH1* h2m(0x0);
1695 TString format;
1696
1697 while ( ( trigger = static_cast<TObjString*>(nextTrigger()) ) )
1698 {
1699 if (bothSigns)
1700 {
1701 format = "/%s/%s/%s/%s/PtEtaMuPlus:py";
1702 }
1703 else
1704 {
1705 format = "/%s/%s/%s/%s/PtEtaMu:py";
1706 }
1707
1708 h1p = hc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
1709
1710 if (!h1p)
1711 {
1712 continue;
1713 }
1714
1715 AliInfoClass(Form("Will use trigger %s",trigger->String().Data()));
1716
1717 h2p = hc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
1718
1719 if ( bothSigns )
1720 {
1721 format.ReplaceAll("Plus","Minus");
1722 h1m = hc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
1723 h2m = hc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
1724 }
1725 else
1726 {
1727 h2m=h2p;
1728 h1m=h1p;
1729 }
1730
1731 if (h1m && h2m && h1p && h2p)
1732 {
1733 runs[trigger->String().Data()].push_back(RunNumberFromFileName(str->String().Data()));
1734 errruns[trigger->String().Data()].push_back(0.5);
1735
1736 double value = 100-h2m->Integral()*100/h1m->Integral();
1737 yminus[trigger->String().Data()].push_back(value);
1738 double e1 = 1.0/TMath::Sqrt(h1m->GetEntries());
1739 double e2 = 1.0/TMath::Sqrt(h2m->GetEntries());
1740 erryminus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
1741
1742 value = 100-h2p->Integral()*100/h1p->Integral();
1743 yplus[trigger->String().Data()].push_back(value);
1744 e1 = 1.0/TMath::Sqrt(h1p->GetEntries());
1745 e2 = 1.0/TMath::Sqrt(h2p->GetEntries());
1746 erryplus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
1747 }
1748 else
1749 {
1750 std::cout << Form("Error : h1m %p h2m %p h1p %p h2p %p",h1m,h2m,h1p,h2p) << std::endl;
1751 }
1752 }
1753
1754 delete hc;
1755 delete cc;
1756 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(str->String().Data()));
1757 delete f;
1758 }
1759
1760 delete triggers;
1761
1762 TFile* f = new TFile(outputFile,outputMode);
1763
1764 std::map<std::string, std::vector<float> >::const_iterator it;
1765
1766 for ( it = runs.begin(); it != runs.end(); ++it )
1767 {
1768 std::string triggerName = it->first;
1769
1770 TGraphErrors* gp = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yplus[triggerName][0],&errruns[triggerName][0],&erryplus[triggerName][0]);
1771 TGraphErrors* gm(0x0);
1772
1773 if ( bothSigns )
1774 {
1775 gm = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yminus[triggerName][0],&errruns[triggerName][0],&erryminus[triggerName][0]);
1776 }
1777
1778 if ( bothSigns )
1779 {
1780 gp->Write(Form("muplus_%s",triggerName.c_str()),TObject::kOverwrite);
1781 gm->Write(Form("muminus_%s",triggerName.c_str()),TObject::kOverwrite);
1782 }
1783 else
1784 {
1785 gp->Write(Form("mu_%s",triggerName.c_str()),TObject::kOverwrite);
1786 }
1787
1788 }
1789
1790 delete f;
1791
1792 return a;
1793}
1794
2754fd07 1795//_____________________________________________________________________________
1796TMap*
1797AliAnalysisMuMu::ComputeJpsiEvolution(const char* filelist, const char* triggerList,
1798 const char* outputFile, const char* outputMode,
1799 Bool_t simulation)
1800{
1801 /// Compute some jpsi information for a list of files / trigger combinations
1802
1803 TObjArray* files = ReadFileList(filelist);
1804
1805 if (!files || files->IsEmpty() ) return 0x0;
1806
1807 TMap results; // one TObjString->TObjArray per file
1808 results.SetOwnerKeyValue(kTRUE,kTRUE);
1809
1810 TIter nextFile(files);
1811 TObjString* str;
1812 UInt_t fitType(0);
1813
1814 while ( ( str = static_cast<TObjString*>(nextFile()) ) )
1815 {
1816 std::cout << str->String().Data() << std::endl;
1817
1818 AliAnalysisMuMu m(str->String().Data());
1819
1820 m.SetDimuonTriggerList(triggerList);
1821 m.SetEventSelectionList("ALL");
1822
1823 TObjArray* array = m.Jpsi(simulation); // the array will contain results for all triggers in fDimuonTriggers variable
1824
1825 if (!array)
1826 {
1827 AliWarningClass(Form("Got no jpsi for %s",str->String().Data()));
1828 }
1829 else
1830 {
1831 Result* r = static_cast<Result*>(array->First());
1832 if (!r) continue;
1833 fitType = r->FitType();
1834 }
1835
1836 results.Add(new TObjString(str->String()), array);
1837 }
1838
1839 if (!results.GetSize()) return 0x0;
1840
1841 // compute the total over all files
1842
1843 TMap* total = new TMap;
1844 total->SetOwnerKeyValue(kTRUE,kTRUE);
1845
1846 nextFile.Reset();
1847 TObjArray* triggers = TString(triggerList).Tokenize(",");
1848
1849 TIter nextTrigger(triggers);
1850 TObjString* trigger(0x0);
1851
1852 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
1853 {
1854 nextFile.Reset();
1855
1856 Int_t nruns(0);
1857 Int_t n(0);
1858 TList l;
1859 Result* ref(0x0);
1860 AliHistogramCollection* hc(0x0);
1861
1862 while ( ( str = static_cast<TObjString*>(nextFile()) ) )
1863 {
1864 TObjArray* a = static_cast<TObjArray*>(results.GetValue(str->String().Data()));
1865
1866 Result* r(0x0);
1867
1868 if (a)
1869 {
1870 r = static_cast<Result*>(a->FindObject(trigger->String().Data()));
1871
1872 if (r)
1873 {
1874 if (!ref) ref = r;
1875
1876 if ( !hc )
1877 {
1878 AliHistogramCollection* htmp = r->HC();
1879 if (!htmp)
1880 {
1881 continue;
1882 }
1883 hc = static_cast<AliHistogramCollection*>(htmp->Clone(Form("hc%d",0)));
1884 }
1885 else
1886 {
1887 l.Add(r->HC());
1888 }
1889
1890 n += r->NofTriggers();
1891 ++nruns;
1892 }
1893 }
1894 }
1895
1896 hc->Merge(&l);
1897
1898 if (!ref) continue;
1899
1900 Result* sum = new Result(ref->TriggerClass(),ref->EventSelection(),
1901 ref->PairSelection(),ref->CentralitySelection(),
1902 n,hc,1,fitType);
1903
1904 sum->SetNofRuns(nruns);
1905
1906 total->Add(new TObjString(trigger->String().Data()),sum);
1907
1908 }
1909
1910 StdoutToAliInfoClass(total->Print(););
1911
1912 TFile* fout = new TFile(outputFile,outputMode);
1913
1914 results.Write("rbr",TObject::kSingleKey|TObject::kOverwrite);
1915
1916 total->Write("total",TObject::kSingleKey|TObject::kOverwrite);
1917
1918 delete fout;
1919
1920 AliInfoClass(Form("%d files analyzed",files->GetEntries()));
1921
1922 return total;
1923}
1924
2f331ac9 1925//_____________________________________________________________________________
1926Bool_t AliAnalysisMuMu::DecodeFileName(const char* filename,
1927 TString& period,
1928 int& esdpass,
1929 int& aodtrain,
1930 int& runnumber)
1931{
1932 esdpass=aodtrain=runnumber=-1;
1933 period="";
1934
1935 TString sfile(gSystem->BaseName(filename));
1936
1937 if (!sfile.BeginsWith("LHC") && !sfile.BeginsWith("SIM") )
1938 {
1939 std::cerr << Form("filename %s does not start with LHC or SIM",filename) << std::endl;
1940 return kFALSE;
1941 }
1942
1943 int year;
1944 char p;
1945 Bool_t ok(kFALSE);
1946
1947 if ( sfile.BeginsWith("LHC") )
1948 {
1949 if ( sfile.Contains("pass") && sfile.Contains("AODMUON") )
1950 {
1951 int pass;
1952 sscanf(sfile.Data(),"LHC%2d%c_pass%d_AODMUON%03d_%09d",&year,&p,&pass,&aodtrain,&runnumber);
1953 period = TString::Format("LHC%2d%c",year,p);
1954 ok=kTRUE;
1955 }
1956 else if ( sfile.Contains("pass") && sfile.Contains("_muon_") && sfile.Contains("AOD000") )
1957 {
1958 // LHC11c_pass2_muon_AOD000_000152087.saf.root
1959 sscanf(sfile.Data(),"LHC%2d%c_pass%d_muon_AOD000_%09d",&year,&p,&esdpass,&runnumber);
1960 period = TString::Format("LHC%2d%c",year,p);
1961 ok=kTRUE;
1962 }
1963 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AODMUON000") )
1964 {
1965 // LHC12h_muon_calo_AODMUON000_000190112.saf.root
1966 sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AODMUON000_%09d",&year,&p,&runnumber);
1967 period = TString::Format("LHC%2d%c",year,p);
1968 ok=kTRUE;
1969 }
1970 else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AOD000") )
1971 {
1972 // LHC12h_muon_calo_AOD000_000190112.saf.root
1973 sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AOD000_%09d",&year,&p,&runnumber);
1974 period = TString::Format("LHC%2d%c",year,p);
1975 ok=kTRUE;
1976 }
1977 else if ( sfile.Contains("AODMUON" ) )
1978 {
1979 sscanf(sfile.Data(),"LHC%2d%c_AODMUON%03d_%09d",&year,&p,&aodtrain,&runnumber);
1980 period = TString::Format("LHC%2d%c",year,p);
1981 ok=kTRUE;
1982 }
1983 else if ( sfile.Contains("AOD000") )
1984 {
1985 sscanf(sfile.Data(),"LHC%2d%c_muons_AOD000_%09d",&year,&p,&runnumber);
1986 period = TString::Format("LHC%2d%c",year,p);
1987 ok=kTRUE;
1988 }
1989 else if ( sfile.Contains("ESD_OUTER000"))
1990 {
1991 sscanf(sfile.Data(),"LHC%2d%c_cpass1_ESD_OUTER000_%09d",&year,&p,&runnumber);
1992 ok=kTRUE;
1993 }
1994 }
1995 else if ( sfile.BeginsWith("SIM_JPSI3" ) )
1996 {
1997 sscanf(sfile.Data(),"SIM_JPSI3_%09d",&runnumber);
1998 ok = kTRUE;
1999 }
2000 else if ( sfile.BeginsWith("SIM_UPSILON" ) )
2001 {
2002 sscanf(sfile.Data(),"SIM_UPSILON_%09d",&runnumber);
2003 ok = kTRUE;
2004 }
2005 else if ( sfile.BeginsWith("SIM_JPSI" ) )
2006 {
2007 sscanf(sfile.Data(),"SIM_JPSI_LHC%2d%c_%09d",&year,&p,&runnumber);
2008 period = TString::Format("LHC%2d%c",year,p);
2009 ok = kTRUE;
2010 }
2011
2012 if (!ok)
2013 {
2014 std::cerr << Form("Can not decode %s",filename) << std::endl;
2015 return kFALSE;
2016 }
2017
2018 return kTRUE;
2019}
2020
2754fd07 2021//______________________________________________________________________________
2022void AliAnalysisMuMu::DrawFill(Int_t run1, Int_t run2, double ymin, double ymax, const char* label)
2023{
2024 AliDebugClass(1,Form("RUN1 %09d RUN2 %09d YMIN %e YMAX %e %s",
2025 run1,run2,ymin,ymax,label));
2026 TBox* b = new TBox(run1*1.0,ymin,run2*1.0,ymax);
2027 b->SetFillColor(5);
2028 b->Draw();
2029 TText* text = new TText((run1+run2)/2.0,ymax*0.85,label);
2030 text->SetTextSize(0.025);
2031 text->SetTextFont(42);
2032 text->SetTextAlign(23);
2033 text->SetTextAngle(45);
2034 text->Draw();
2035}
2036
2f331ac9 2037//_____________________________________________________________________________
2038TString
2039AliAnalysisMuMu::ExpandPathName(const char* file)
2040{
2041 // An expand method that lives alien URL as they are
2042 TString sfile;
2043
2044 if ( !sfile.BeginsWith("alien://") )
2045 {
2046 return gSystem->ExpandPathName(file);
2047 }
2048 else
2049 {
2050 if (!gGrid) TGrid::Connect("alien://");
2051 if (!gGrid) return "";
2052 }
2053
2054 return file;
2055}
2056
2057//_____________________________________________________________________________
2058TFile*
2059AliAnalysisMuMu::FileOpen(const char* file)
2060{
2061 // Open a file after expansion of its name
2062
2063 return TFile::Open(ExpandPathName(file).Data());
2064}
2065
2066//_____________________________________________________________________________
2067ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
2068{
2069 // Get the expected (from OCDB scalers) trigger count
2070
2071 AliAnalysisTriggerScalers ts(runNumber,fgOCDBPath.Data());
2072
2754fd07 2073 TObjArray* triggers = TString(triggerList).Tokenize(",");
2f331ac9 2074 TObjString* trigger;
2075 TIter next(triggers);
2076 ULong64_t n(0);
2077
2078 while ( ( trigger = static_cast<TObjString*>(next()) ) )
2079 {
2080 AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
2081 if (item)
2082 {
2083 n += item->Value();
2084 }
2085 delete item;
2086 }
2087 delete triggers;
2088
2089 return n;
2090}
2091
2092//_____________________________________________________________________________
2093UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList, const char* eventSelection, Int_t runNumber)
2094{
2095 TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
2096 TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
2097 TObjArray* a = TString(triggerList).Tokenize(" ");
2098 TIter next(a);
2099 TObjString* str;
2100
2101 UInt_t n(0);
2102
2103 TString sEventSelection(eventSelection);
2104 sEventSelection.ToUpper();
2105
2106 if ( kevent->FindObject(sEventSelection.Data()) )
2107 {
2108 while ( ( str = static_cast<TObjString*>(next()) ) )
2109 {
2110 if ( ktrigger->FindObject(str->String().Data()) )
2111 {
2112 if ( runNumber < 0 )
2113 {
2114 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
2115 }
2116 else
2117 {
2118 n += static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
2119 }
2120 }
2121 }
2122 }
2123
2124 delete a;
2125 delete ktrigger;
2126 delete kevent;
2127 return n;
2128}
2129
2130//_____________________________________________________________________________
2131TObjArray*
2132AliAnalysisMuMu::Jpsi(Bool_t simulation)
2133{
2754fd07 2134 // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
2135
2f331ac9 2136 TObjArray* a = new TObjArray;
2137 a->SetOwner(kTRUE);
2138
2139 if ( simulation )
2140 {
2141 const char* selection = "ALL";
2142
2754fd07 2143 Add(a,GetResult(*fHistogramCollection,*fCounterCollection,"ANY",selection,"pMATCHLOWRABSBOTH","PP",Result::kJpsi));
2f331ac9 2144 }
2145 else
2146 {
2754fd07 2147 TObjArray* triggerArray = fDimuonTriggers.Tokenize(",");
2148 TObjArray* eventTypeArray = fEventSelectionList.Tokenize(",");
2149 TObjArray* pairCutArray = fPairSelectionList.Tokenize(",");
2f331ac9 2150
2754fd07 2151 TIter nextTrigger(triggerArray);
2152 TIter nextEventType(eventTypeArray);
2153 TIter nextPairCut(pairCutArray);
2154
2155 TObjString* trigger;
2156 TObjString* eventType;
2157 TObjString* pairCut;
2158
2159 while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
2f331ac9 2160 {
2754fd07 2161 AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
2162
2163 nextEventType.Reset();
2164
2165 while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
2f331ac9 2166 {
2754fd07 2167 AliDebug(1,Form("EVENTTYPE %s",eventType->String().Data()));
2f331ac9 2168
2754fd07 2169 nextPairCut.Reset();
2f331ac9 2170
2754fd07 2171 while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
2172 {
2173 AliDebug(1,Form("PAIRCUT %s",pairCut->String().Data()));
2174 Add(a,GetResult(*fHistogramCollection,*fCounterCollection,
2175 trigger->String().Data(),
2176 eventType->String().Data(),
2177 pairCut->String().Data(),
2178 "PP",
2179 Result::kJpsi | Result::kJpsiPrime,2));
2180 }
2f331ac9 2181 }
2182 }
2183
2754fd07 2184 delete triggerArray;
2185 delete eventTypeArray;
2186 delete pairCutArray;
2f331ac9 2187 }
2188
2189 if ( a->GetLast() < 0 )
2190 {
2191 delete a;
2192 a = 0;
2193 }
2194 return a;
2195
2196}
2197
2198//_____________________________________________________________________________
2754fd07 2199void
2200AliAnalysisMuMu::PlotJpsiEvolution(const char* resultFile, const char* triggerList, Bool_t fillBoundaries,
2201 const char* efficiencyFile)
2202{
2203 if ( efficiencyFile && strlen(efficiencyFile) > 0 )
2204 {
2205 std::ifstream in(gSystem->ExpandPathName(efficiencyFile));
2206 if (!in.bad())
2207 {
2208 char line[1024];
2209 int run;
2210 float eff, error;
2211 while ( in.getline(line,1023,'\n') )
2212 {
2213 sscanf(line,"%d %f ± %f",&run,&eff,&error);
2214 AliInfoClass(Form("%09d %8.6f ± %8.6f",run,eff,error));
2215 }
2216 }
2217
2218 return;
2219 }
2f331ac9 2220
2754fd07 2221 TFile* f = TFile::Open(gSystem->ExpandPathName(resultFile));
2f331ac9 2222
2754fd07 2223 std::map<int, std::pair<int,int> > fills;
2f331ac9 2224
2754fd07 2225 TMap* m = static_cast<TMap*>(f->Get("rbr"));
2226
2227 TIter next(m);
2f331ac9 2228 TObjString* str;
2f331ac9 2229
2754fd07 2230 TObjArray files;
2231 files.SetOwner(kTRUE);
2232
2233 while ( ( str = static_cast<TObjString*>(next())) )
2f331ac9 2234 {
2754fd07 2235 files.Add(new TObjString(str->String()));
2f331ac9 2236 }
2237
2754fd07 2238 files.Sort();
2f331ac9 2239
2754fd07 2240 std::map<std::string, std::vector<float> > x_jpsirate;
2241 std::map<std::string, std::vector<float> > y_jpsirate;
2242 std::map<std::string, std::vector<float> > xerr_jpsirate;
2243 std::map<std::string, std::vector<float> > yerr_jpsirate;
2f331ac9 2244
2754fd07 2245 TIter nextTrigger(TString(triggerList).Tokenize(","));
2246 TObjString* trigger(0x0);
2f331ac9 2247
2754fd07 2248 int runMin(100000000);
2249 int runMax(0);
2250
2251 TIter nextFile(&files);
2252
2253 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2f331ac9 2254 {
2754fd07 2255 TString triggerClass(trigger->String());
2f331ac9 2256
2754fd07 2257 nextFile.Reset();
2f331ac9 2258
2f331ac9 2259
2754fd07 2260 while ( ( str = static_cast<TObjString*>(nextFile())) )
2f331ac9 2261 {
2754fd07 2262 TObjArray* a = static_cast<TObjArray*>(m->GetValue(str->String().Data()));
2263 if (!a) continue;
2264 Result* r = static_cast<Result*>(a->FindObject(triggerClass.Data()));
2265 if (!r) continue;
2266
2267 TString period;
2268 int aodtrain,esdpass,runnumber;
2269
2270 if ( DecodeFileName(str->String().Data(),period,esdpass,aodtrain,runnumber) )
2f331ac9 2271 {
2754fd07 2272 runMin = TMath::Min(runMin,runnumber);
2273 runMax = TMath::Max(runMax,runnumber);
2f331ac9 2274
2754fd07 2275 x_jpsirate[triggerClass.Data()].push_back(runnumber);
2276 xerr_jpsirate[triggerClass.Data()].push_back(0.5);
2277
2278 if ( fillBoundaries )
2f331ac9 2279 {
2754fd07 2280 AliAnalysisTriggerScalers ts(runnumber,fgOCDBPath.Data());
2281 int fill = ts.GetFillNumberFromRunNumber(runnumber);
2282
2283 if (fills.count(fill))
2f331ac9 2284 {
2754fd07 2285 std::pair<int,int>& p = fills[fill];
2286 p.first = TMath::Min(runnumber,p.first);
2287 p.second = TMath::Max(runnumber,p.second);
2f331ac9 2288 }
2289 else
2290 {
0775d258 2291 fills[fill] = std::make_pair<int,int>(runnumber,runnumber);
2f331ac9 2292 }
2293 }
2f331ac9 2294
2754fd07 2295 Double_t y(0.0);
2296 Double_t yerr(0.0);
2f331ac9 2297
2754fd07 2298 if ( TMath::Finite(r->SigmaJpsi()) && r->NofTriggers() > 10 )
2f331ac9 2299 {
2754fd07 2300 y = 100*r->NofJpsi()/r->NofTriggers();
2f331ac9 2301
2754fd07 2302 if ( r->NofJpsi() > 0 )
2303 {
2304 yerr = y * TMath::Sqrt( (r->ErrorOnNofJpsi()*r->ErrorOnNofJpsi())/(r->NofJpsi()*r->NofJpsi()) + 1.0/r->NofTriggers());
2305 }
2f331ac9 2306 }
2754fd07 2307
2308 y_jpsirate[triggerClass.Data()].push_back(y);
2309 yerr_jpsirate[triggerClass.Data()].push_back(yerr);
2f331ac9 2310 }
2f331ac9 2311 }
2f331ac9 2312 }
2f331ac9 2313
2754fd07 2314 delete f;
2f331ac9 2315
2754fd07 2316 TCanvas* c = new TCanvas("cJpsiRateEvolution","cJpsiRateEvolution");
2f331ac9 2317
2754fd07 2318 c->Draw();
2f331ac9 2319
2754fd07 2320 Double_t ymin(0);
2321 Double_t ymax(2);
2f331ac9 2322
2754fd07 2323 TH2* h = new TH2F("h","h;RunNumber;J/#psi per CMUL (%)",100,runMin,runMax,100,ymin,ymax);
2f331ac9 2324
2754fd07 2325 gStyle->SetOptTitle(0);
2326 gStyle->SetOptStat(0);
2f331ac9 2327
2754fd07 2328 h->GetXaxis()->SetNoExponent();
2f331ac9 2329
2754fd07 2330 h->Draw();
2331
2332 if (fillBoundaries)
2333 {
2334 std::map<int, std::pair<int,int> >::const_iterator it;
2335
2336 for ( it = fills.begin(); it != fills.end(); ++it )
2337 {
2338 const std::pair<int,int>& p = it->second;
2339 TString fillnumber;
2340 fillnumber.Form("%d",it->first);
2341 DrawFill(p.first,p.second,ymin,ymax,fillnumber.Data());
2342 }
2343 }
2344
2345 h->Draw("sameaxis");
2f331ac9 2346
2754fd07 2347 //c->RedrawAxis("g");
2348
2349 nextTrigger.Reset();
2f331ac9 2350
2754fd07 2351 int i(0);
2352 int color[] = { 1,2,4,5,6 };
2353 int marker[] = { 20,23,25,21,22 };
2f331ac9 2354
2754fd07 2355 while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2356 {
2357 std::vector<float>& x = x_jpsirate[trigger->String().Data()];
2358 std::vector<float>& y = y_jpsirate[trigger->String().Data()];
2359 std::vector<float>& xerr = xerr_jpsirate[trigger->String().Data()];
2360 std::vector<float>& yerr = yerr_jpsirate[trigger->String().Data()];
2361
2362 TGraphErrors* g = new TGraphErrors(x.size(),&x[0],&y[0],&xerr[0],&yerr[0]);
2363
2364 g->SetLineColor(color[i]);
2365 g->SetMarkerColor(color[i]);
2366 g->SetMarkerStyle(marker[i]);
2367 g->GetXaxis()->SetNoExponent();
2368 g->Draw("LP");
2369// g->Print();
2370
2371 Double_t m2 = g->GetMean(2);
2372
2373 TLine* line = new TLine(runMin,m2,runMax,m2);
2374 line->SetLineColor(color[i]);
2375 line->Draw();
2376
2377 AliInfoClass(Form("TRIGGER %s MEAN %7.2f",trigger->String().Data(),m2));
2378 ++i;
2379 }
2f331ac9 2380
2f331ac9 2381
2f331ac9 2382}
2383
2384//_____________________________________________________________________________
2754fd07 2385AliAnalysisMuMu::Result*
2f331ac9 2386AliAnalysisMuMu::GetResult(const AliHistogramCollection& hc,
2754fd07 2387 AliCounterCollection& cc,
2388 const char* base,
2389 const char* selection,
2390 const char* paircut,
2391 const char* centrality,
2392 UInt_t fitType,
2393 Int_t nrebin)
2f331ac9 2394{
2395 Result* r(0x0);
2754fd07 2396
2f331ac9 2397 TString trigger = FindTrigger(hc,base,selection,paircut,centrality);
2398
2754fd07 2399 if ( trigger == "" )
2400 {
2401 return 0;
2402 }
2f331ac9 2403
2404 Int_t ntrigger = (Int_t)cc.GetSum(Form("trigger:%s/event:%s",trigger.Data(),selection));
2405
2406// new TCanvas;
2407
2754fd07 2408 r = new Result(trigger.Data(),
2409 selection,
2410 paircut,
2411 centrality,
2f331ac9 2412 ntrigger,
2754fd07 2413 hc.Project(selection,trigger,centrality,paircut),
2f331ac9 2414 nrebin,
2415 fitType);
2416
2417 return r;
2418}
2419
2420//_____________________________________________________________________________
2421Bool_t
2422AliAnalysisMuMu::GetCollections(const char* rootfile,
2423 AliHistogramCollection*& hc,
2424 AliCounterCollection*& cc)
2425{
2426 hc = 0x0;
2427 cc = 0x0;
2428
2429 TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
2430
2431 if (!f)
2432 {
2433 f = TFile::Open(rootfile);
2434 }
2435
2436 if ( !f || !f->IsOpen() )
2437 {
2438 return kFALSE;
2439 }
2440
2441 TList* list = static_cast<TList*>(f->Get("chist"));
2442
2443 if (!list) return kFALSE;
2444
2445 hc = static_cast<AliHistogramCollection*>(list->At(0));
2446 cc = static_cast<AliCounterCollection*>(list->At(1));
2447
2448 return kTRUE;
2449}
2450
2451//_____________________________________________________________________________
2452void AliAnalysisMuMu::PlotBackgroundEvolution(const char* gfile, const char* triggerList)
2453{
2454 // plot the graphs found in the file (and generated using the ComputeBackgroundEvolution() method)
2455
2456 TFile* f = TFile::Open(ExpandPathName(gfile).Data());
2457
2458 if ( !f || !f->IsOpen() )
2459 {
2460 return;
2461 }
2462
2463 SetColorScheme();
2464
2465
2466 TCanvas* c = new TCanvas("background-evolution","background-evolution");
2467
2468 c->Draw();
2469
2470 TLegend* l = new TLegend(0.4,0.6,0.97,0.97);
2471 l->SetFillColor(0);
2472 l->SetTextColor(AliAnalysisMuMu::kBlue);
2473 l->SetLineColor(AliAnalysisMuMu::kBlue);
2474
2754fd07 2475 TObjArray* triggers = TString(triggerList).Tokenize(",");
2f331ac9 2476
2477 gStyle->SetOptTitle(0);
2478
2479 TObjString* str(0x0);
2480 TIter next(triggers);
2481 Int_t i(0);
2482 Int_t run1(99999999);
2483 Int_t run2(0);
2484
2485 while ( ( str = static_cast<TObjString*>(next()) ) )
2486 {
2487 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2488 if (!g) continue;
2489 run1 = TMath::Min(run1,TMath::Nint(g->GetX()[0]));
2490 run2 = TMath::Max(run2,TMath::Nint(g->GetX()[g->GetN()-1]));
2491 }
2492
2493 AliInfoClass(Form("run1 %d run2 %d",run1,run2));
2494
2495 TH2* hframe = new TH2F("hframe","hframe",run2-run1+1,run1,run2,100,0,100);
2496 hframe->Draw();
2497 hframe->GetXaxis()->SetNoExponent();
2498 hframe->GetYaxis()->SetTitle("Background percentage");
2499
2500 next.Reset();
2501
2502 while ( ( str = static_cast<TObjString*>(next()) ) )
2503 {
2504 TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2505 if (!g)
2506 {
2507 AliErrorClass(Form("Graph mu_%s not found",str->String().Data()));
2508 continue;
2509 }
2510
2511 Int_t color(i+1);
2512
2513 if (i==0) color = AliAnalysisMuMu::kBlue;
2514 if (i==1) color = AliAnalysisMuMu::kOrange;
2515
2516 g->SetLineColor(color);
2517 g->SetMarkerColor(color);
2518 g->SetMarkerStyle(20+i);
2519
2520 g->Draw("LP");
2521
2522 TLegendEntry* le = l->AddEntry(g,str->String().Data(),"lp");
2523 le->SetTextColor(color);
2524
2525 g->GetYaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2526 g->GetXaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2527// g->GetXaxis()->SetNoExponent();
2528// g->GetYaxis()->SetRangeUser(0,100);
2529// g->GetYaxis()->SetTitle("Background percentage");
2530 g->Print();
2531
2532 ++i;
2533 }
2534
2535 l->Draw();
09d5920f 2536 delete triggers;
2f331ac9 2537}
2538
2539//_____________________________________________________________________________
2540void AliAnalysisMuMu::SetColorScheme()
2541{
2542 new TColor(AliAnalysisMuMu::kBlue,4/255.0,44/255.0,87/255.0,"my blue");
2543 new TColor(AliAnalysisMuMu::kOrange,255/255.0,83/255.0,8/255.0,"my orange");
2544 new TColor(AliAnalysisMuMu::kGreen,152/255.0,202/255.0,52/255.0,"my green");
2545
2546 gStyle->SetGridColor(AliAnalysisMuMu::kBlue);
2547
2548 gStyle->SetFrameLineColor(AliAnalysisMuMu::kBlue);
2549 gStyle->SetAxisColor(AliAnalysisMuMu::kBlue,"xyz");
2550 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue,"xyz");
2551
2552 gStyle->SetTitleColor(AliAnalysisMuMu::kBlue);
2553 gStyle->SetTitleTextColor(AliAnalysisMuMu::kBlue);
2554 gStyle->SetLabelColor(AliAnalysisMuMu::kBlue);
2555 gStyle->SetStatTextColor(AliAnalysisMuMu::kBlue);
2556
2557 gStyle->SetOptStat(0);
2558}
2559
2560//_____________________________________________________________________________
2561void
2562AliAnalysisMuMu::SinglePtPlot(const char* rootfile)
2563{
2564 AliHistogramCollection* histogramCollection(0x0);
2565 AliCounterCollection* counterCollection(0x0);
2566
2567 if (!GetCollections(rootfile,histogramCollection,counterCollection))
2568 {
2569 return;
2570 }
2571
2572 TCanvas* c1 = new TCanvas("singlept","singlept");
2573
2574 gStyle->SetTextSize(1.0);
2575
2576 c1->SetLogy();
2577
2578 TLegend* l = new TLegend(0.12,0.12,0.4,0.3,"","NDC");
2579 l->SetFillStyle(0);
2580 l->SetLineWidth(0);
2581 l->SetLineColor(0);
2582 l->SetTextColor(kBlack);
2583 l->SetTextFont(42);
2584 l->SetTextSize(0.035);
2585
2586 const char* cuts[] = { "ALL","ETA","ETARABS","ETARABSMATCH","ETARABSMATCHDCA" };
2587 const char* cutnames[] = {
2588 "all",
2589 "+ -4 < #eta < -2.5",
2590 "+ 171^{#circ} < #theta_{abs} < 178^{#circ}",
2591 "+ trigger matching",
2592 "+ PxDCA"
2593 };
2594 const int colors[] = { 1,2,3,4,6 };
2595
2596 for ( int i = 0; i < 5; ++i )
2597 {
2598 TH1* minus = histogramCollection->Histo(Form("/PS/CPBI1MSH-B-NOPF-MUON/CENT80/s%s/PtEtaMuMinus:py",cuts[i]));
2599 TH1* plus = histogramCollection->Histo(Form("/PS/CPBI1MSH-B-NOPF-MUON/CENT80/s%s/PtEtaMuPlus:py",cuts[i]));
2600 if (!minus || !plus)
2601 {
2602 AliErrorClass(Form("Form cannot get histos for cut %s",cuts[i]));
2603 continue;
2604 }
2605 TH1* h = static_cast<TH1*>(minus->Clone(Form("h%d",i)));
2606 h->Add(plus);
2607 h->SetDirectory(0);
2608 h->SetLineColor(colors[i]);
2609 h->SetMinimum(1);
2610 h->SetMarkerSize(0);
2611 h->SetLineWidth(2);
2612 h->SetStats(0);
2613 h->SetBit(TH1::kNoTitle);
2614 if (i==0)
2615 {
2616 h->Draw("e");
2617 }
2618 else
2619 {
2620 h->Draw("esame");
2621 }
2622 h->GetYaxis()->SetTitle("dN/dp_{t} (counts/0.5 GeV/c)");
2623 h->GetXaxis()->SetTitle("p_{t} (GeV/c)");
2624
2625 h->GetXaxis()->SetLabelSize(0.04);
2626 h->GetYaxis()->SetLabelSize(0.04);
2627 h->GetYaxis()->SetTitleSize(0.04);
2628 h->GetXaxis()->SetTitleSize(0.04);
2629
2630 h->GetXaxis()->SetRangeUser(0,30);
2631 l->AddEntry(h,cutnames[i]);
2632 }
2633
2634 l->Draw();
2635
2636 TList extralines;
2637 extralines.SetOwner(kTRUE);
2638
2639 extralines.Add(new TObjString("PbPb #sqrt{s_{NN}}=2.76 TeV"));
2640 extralines.Add(new TObjString("MUON high p_{t} trigger events"));
2641 extralines.Add(new TObjString("w/ phys. sel."));
2642 extralines.Add(new TObjString("w/ reco. vertex"));
2643 extralines.Add(new TObjString("centrality 0-80 %"));
2644
2645 ALICEseal(0.5,0.93,&extralines);
2646
2647 delete histogramCollection;
2648 delete counterCollection;
2649}
2650
2651//_____________________________________________________________________________
2754fd07 2652void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList, Bool_t compact)
2f331ac9 2653{
2654 // Give the fraction of triggers (in triggerList) relative
2655 // to what is expected in the scalers
2656
2657 TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
2658
2659 AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
2660
2661 AliLog::SetGlobalLogLevel(AliLog::kFatal);
2662
2663 if (!fHistogramCollection || !fCounterCollection) return;
2664
2665 TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2666 TIter nextRun(runs);
2667
2668 TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2669 TIter nextTrigger(triggers);
2670
2671 TObjString* srun;
2672 TObjString* strigger;
2673
2674 TString striggerList(triggerList);
2675
2676 ULong64_t total(0);
2677 ULong64_t totalExpected(0);
2678
2679 while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2680 {
2681 std::cout << Form("RUN %09d ",srun->String().Atoi());
2682
2754fd07 2683 if (!compact) std::cout << std::endl;
2684
2f331ac9 2685 nextTrigger.Reset();
2686
2687 while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2688 {
2689 if ( !striggerList.Contains(strigger->String().Data()) )
2690 {
2691 continue;
2692 }
2693 ULong64_t n = fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2694 strigger->String().Data(),"ALL",srun->String().Atoi()));
2695
2696 ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
2697
2698
2699 total += n;
2700 totalExpected += expected;
2701
2754fd07 2702 std::cout << Form("%30s %9lld expected %9lld ",strigger->String().Data(),n,expected);
2f331ac9 2703
2704 if ( expected > 0 ) {
2705 std::cout << Form("fraction %5.1f %%",n*100.0/expected);
2706 }
2754fd07 2707 if (!compact)
2708 {
2709 std::cout << std::endl;
2710 }
2f331ac9 2711 }
2712 std::cout << std::endl;
2713 }
2714
2715 std::cout << Form("TOTAL %lld expected %lld fraction %5.1f %%",
2716 total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
2717
2718 AliLog::SetGlobalLogLevel(oldLevel);
09d5920f 2719 delete triggers;
2720 delete runs;
2f331ac9 2721}
2722
2723//_____________________________________________________________________________
2724void AnalyisResultLocation(const char* runlist, const char* basedir, const char* what)
2725{
2726 std::ifstream in(runlist);
2727 int run;
2728 while ( in >> run )
2729 {
2730 std::cout << Form("%s/%09d/%s",basedir,run,what) << std::endl;
2731 }
2732}
2733
2734//_____________________________________________________________________________
2735void plot(const char* file="results.root", const char* pdfname="toto.pdf")
2736{
2737 gROOT->SetStyle("Plain");
2738 gStyle->SetOptTitle(0);
2739
2740 TFile* f = TFile::Open(file);
2741
2742 TObjArray* a = static_cast<TObjArray*>(f->Get("results"));
2743
2744 if (!a) return;
2745
2746 TCanvas* c = new TCanvas("jpsiresults","jpsiresults");
2747
2748 c->Draw();
2749
2750 TLegend* l = new TLegend(0.5,0.7,0.9,0.9);
2751 l->SetFillColor(0);
2752 TH1* h0(0x0);
2753
2754 for ( int i = 0; i <= a->GetLast(); ++i )
2755 {
2756 AliAnalysisMuMu::Result* r = static_cast<AliAnalysisMuMu::Result*>(a->At(i));
2757 r->Print();
2758
2759 TH1* h = static_cast<TH1*>(r->Minv());
2760
2761 h->SetStats(kFALSE);
2762 h->SetXTitle("M_{#mu^{+}#mu^{-}} (Gev/c^{2})");
2763 h->SetLineColor(i+1);
2764
2765 h->SetMaximum(5E3);
2766 h->SetMarkerStyle(0);
2767
2768 if ( i == 0 )
2769 {
2770 h0 = h;
2771 h->Draw("hist");
2772 }
2773 else
2774 {
2775 h->Draw("histsame");
2776 }
2777
2778 TObjArray* n = TString(r->TriggerClass()).Tokenize("-");
2779 TObjString* nn = static_cast<TObjString*>(n->First());
2780 l->AddEntry(h,Form("%6d %s",r->NofTriggers(),nn->String().Data()));
2781
2782 delete n;
2783 }
2784
2785 h0->Draw("histsame");
2786 l->Draw();
2787 c->SetLogy();
2788
2789 delete f;
2790
2791 c->SaveAs(pdfname);
2792}
2793
2794//_____________________________________________________________________________
2795TObjArray*
2796AliAnalysisMuMu::ReadFileList(const char* filelist)
2797{
2798 //
2799 // read the filelist and try to order it by runnumber
2800 //
2801 // filelist can either be a real filelist (i.e. a text file containing
2802 // root filenames) or a root file itself.
2803 //
2804
2805 char line[1024];
2806
2807 TObjArray* files = new TObjArray;
2808 files->SetOwner(kTRUE);
2809
2810 TString sfilelist(ExpandPathName(filelist));
2811
2812 if ( sfilelist.EndsWith(".root") )
2813 {
2814 files->Add(new TObjString(sfilelist.Data()));
2815 return files;
2816 }
2817
2818 std::set<int> runnumbers;
2819 std::map<int,std::string> filemap;
2820
2821 std::ifstream in(sfilelist.Data());
2822
2823 TString period;
2824 int aodtrain,esdpass,runnumber;
2825
2826 while ( in.getline(line,1022,'\n') )
2827 {
2828 DecodeFileName(line,period,esdpass,aodtrain,runnumber);
2829
2830 AliDebugClass(1,Form("line %s => period %s esdpass %d aodtrain %d runnumber %09d",
2831 line,period.Data(),esdpass,aodtrain,runnumber));
2832
2833 filemap.insert(std::make_pair<int,std::string>(runnumber,line));
2834 runnumbers.insert(runnumber);
2835 }
2836
2837 in.close();
2838
2839 std::set<int>::const_iterator it;
2840
2841 for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2842 {
2843 files->Add(new TObjString(filemap[*it].c_str()));
2844 }
2845
2846 return files;
2847}
2848
2849//_____________________________________________________________________________
2850Int_t
2851AliAnalysisMuMu::RunNumberFromFileName(const char* filename)
2852{
2853 TString period;
2854 int esdpass,aodtrain,runnumber;
2855 Bool_t ok = DecodeFileName(filename,period,esdpass,aodtrain,runnumber);
2856 if ( ok ) return runnumber;
2857 return -1;
2858}