Adding PiKP-only histograms and eliminating a number of switches where histograms...
[u/mrichter/AliRoot.git] / MUON / MUONplotefficiency.C
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 /* $Id$ */
17
18 /// \ingroup macros
19 /// \file MUONplotefficiency.C
20 /// \brief Macro to compare/plot the Ntuple stored in MUONefficiency.root
21 ///
22 /// Comparison is done between the generated and reconstructed resonance.
23 /// Default is Upsilon but Jpsi can be studied if the ResType argument is changed.
24 /// This allows to determine several important quantities: 
25 /// - reconstruction efficiency (vs pt,y), invariant mass peak variations (vs pt,y)
26 /// - reconstructed tracks and trigger tracks  matching efficiency
27 ///
28 /// \author Christophe Suire, IPN Orsay
29
30 #if !defined(__CINT__) || defined(__MAKECINT__)
31 // ROOT includes
32 #include "TStyle.h"
33 #include "TNtuple.h"
34 #include "TF1.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TLatex.h"
38 #include "TFile.h"
39 #include "TLegend.h"
40 #include "TCanvas.h"
41 #include "TGraphErrors.h"
42 #include <Riostream.h>
43 #include "TMath.h"
44
45 #endif
46
47
48
49 Double_t fitlandau_a(Double_t *x, Double_t *par)
50 {
51   Float_t val = -(x[0] - TMath::Abs(par[1])) / TMath::Abs(par[2]);
52   Double_t result = 0.399 *TMath::Abs(par[0]) * exp( -0.5* (val + exp(-val)) );
53   return result;
54 }
55
56 Double_t fitlandau(Double_t *x, Double_t *par)
57 {
58   Float_t val = -x[0] + 2*par[1]; // - TMath::Abs(par[1]) / TMath::Abs(par[2]);
59   Double_t result = par[0] * TMath::Landau(val,par[1],par[2]);
60   return result;
61 }
62 Double_t fitlandaugauss(Double_t *x,Double_t *par)
63 {
64   //cout << "x[0] : " << x[0] << endl ;
65
66    Float_t min=par[4];  Float_t max=par[5]; 
67
68   Int_t NStep;
69   Float_t Step=.1;
70   NStep=(Int_t)((max-min)/Step)+1;
71   Float_t tx;
72   Double_t result=0;
73   for(Int_t iStep=0;iStep<NStep;iStep++){
74     tx=min+iStep*Step;
75     result+=par[0]*TMath::Landau(-tx+2*par[1],par[1],par[2])
76       *TMath::Gaus(x[0]-tx,0.,par[3])
77       *Step;
78     //cout <<"x[0]-tx/2 : " << x[0] - tx << endl ;
79   }
80   return result;
81 }
82
83
84 Int_t MUONplotefficiency(Int_t ResType = 553, Int_t fittype = 1){
85  
86   Int_t SAVE=1;
87   Int_t WRITE=1;
88   
89   //Double_t MUON_MASS = 0.105658369;
90   Double_t UPSILON_MASS = 9.4603 ;
91   Double_t JPSI_MASS = 3.096916;
92   //Double_t PI = 3.14159265358979312; 
93   Double_t RESONANCE_MASS = 0.; 
94
95   if (ResType==553)  
96      RESONANCE_MASS = UPSILON_MASS;
97   if (ResType==443)  
98      RESONANCE_MASS = JPSI_MASS ;
99     
100   // 553 is the pdg code for Upsilon 
101   // 443 is the pdg code for Jpsi 
102
103   Int_t fitfunc = fittype;
104   // 0 gaussian fit +/- 450 MeV/c2
105   // 1 gaussian fit +/- 250 MeV/c2
106   // 2 approximated landau fit reverted 
107   // 3 landau fit reverted convoluated with a gaussian 
108   // 4 reverted landau fitlan
109
110   
111   Float_t gaussianFitWidth = 0.450 ; // in GeV/c2
112   if (fitfunc==1) gaussianFitWidth = 0.250 ; 
113   
114   // fit range : very important for LandauGauss Fit
115   Float_t FitRange = 0.9 ; //GeV/c2
116   Float_t FitLow = 0.0  ;  Float_t FitHigh = 100. ;
117   FitLow  = RESONANCE_MASS-FitRange ;
118   FitHigh = RESONANCE_MASS+FitRange ;
119
120   
121   //gStyle->Reset();
122   // two ways to set the stat box for all histos
123   gStyle->SetOptStat(1110);
124   
125   // place the stat box
126   gStyle->SetStatX(0.5);
127   gStyle->SetStatY(0.550);
128   //gStyle->SetStatW(0.2);
129   gStyle->SetStatH(0.2);
130   
131   // choose histos with stat
132   // gStyle->SetOptStat(1111);
133   //hist->SetStats(kTRUE); or hist->SetStats(kFALSE);
134   
135   //gStyle->SetOptFit(1111);
136   gStyle->SetOptFit(1);
137   //gStyle->SetOptFit(0);
138
139   gStyle->SetOptTitle(0);
140   gStyle->SetPalette(1,0);
141   //gStyle->SetOptLogy(1);
142   //gStyle->SetOptLogx(0);
143   
144   gStyle->SetFrameFillColor(10);
145   gStyle->SetFillColor(10);
146   gStyle->SetCanvasColor(10);
147
148
149   TFile *effFile = new TFile("MUONefficiency.root","READ");
150   TNtuple *Ktuple = (TNtuple*)effFile->Get("Ktuple");
151   //("Ktuple","Kinematics NTuple","ev:npart:id:idmo:idgdmo:p:pt:y:theta:pseudorap:vx:vy:vz");
152   TNtuple *ESDtuple = (TNtuple*)effFile->Get("ESDtuple");
153   //("ESDtuple","Reconstructed Mu+Mu- pairs","ev:pt:y:theta:minv:pt1:y1:theta1:q1:pt2:y2:theta2:q2");
154   TNtuple *ESDtupleBck = (TNtuple*)effFile->Get("ESDtupleBck");
155   //("ESDtuple","Reconstructed Mu+Mu- pairs","ev:pt:y:theta:minv:pt1:y1:theta1:q1:pt2:y2:theta2:q2");
156
157
158   /*********************************/
159   // Histograms limits and binning
160   Float_t ptmin = 0.0;      Float_t ptmax =  20.0;     Int_t ptbins = 10; 
161   Float_t ymin = -4.5;       Float_t ymax =  -2.;       Int_t ybins = 10;
162   // Float_t thetamin = 165.;  Float_t thetamax = 180.;   Int_t thetabins = 100;  
163   
164   // Float_t etacutmin = -4.04813 ;
165   // Float_t etacutmax = -2.54209 ;
166
167   Float_t invMassMin = .0 ;   Float_t invMassMax = 15.0 ;   Int_t   invMassBins = 100 ;
168
169   if (ResType==443){
170     invMassMin = 1.0 ;  invMassMax = 4.5 ;  
171     ptmin = 0.0;     ptmax =  10.0; 
172   }
173   if (ResType==553){
174     invMassMin =  7.0 ;  invMassMax = 11.0 ;  
175     ptmin = 0.0;   ptmax =  20.0; 
176   }
177   
178   /*********************************/
179   // Values used to define the acceptance
180   //
181   // Float_t thetacutmin = 171.0;  
182   // Float_t thetacutmax = 178.;
183   
184   Float_t ptcutmin = 0.;
185   Float_t ptcutmax = 20.;
186
187   if (ResType==443){
188     ptcutmin = 0.; ptcutmax = 10.;
189   }
190
191   if (ResType==553){
192     ptcutmin = 0.; ptcutmax = 20.;
193   }
194   
195
196   Float_t masssigma  = 0.1; // 100 MeV/c2 is a correct estimation 
197   Float_t masscutmin = 0; Float_t masscutmax = 0 ;    
198   if (masssigma){
199     masscutmin = RESONANCE_MASS - 3.0*masssigma ; 
200     masscutmax = RESONANCE_MASS + 3.0*masssigma ;
201   }
202   else {
203     masscutmin = RESONANCE_MASS - 1.0 ; 
204     masscutmax = RESONANCE_MASS + 1.0 ;
205   }
206
207
208   // here no cut on theta is necesary since during the simulation 
209   // gener->SetChildThetaRange(171.0,178.0);
210   // thus the resonance are generated in 4 PI but only the ones for which the decay muons 
211   // are in the dimuon arm acceptance
212   // the pt cut is also "critical", in the generation part, resonance are generated within  the range 
213   // (ptcutmin,ptcutmax). During the reconstruction, the resonance could have a pt greater than ptcutmax !!
214   // Should these resonances be cut or not ?
215   // probably not but they will be for now (~ 1/2000)
216
217   // Another acceptance cut to add is a range for the  recontructed  invariant mass, it is 
218   // obvious that an upsilon reconstructed with a mass  of 5 GeV/c2 is not correct. Thus
219   Char_t ResonanceAccCutMC[200];  
220   sprintf(ResonanceAccCutMC,"pt>= %.2f && pt<= %.2f",ptcutmin,ptcutmax);
221
222   Char_t ResonanceAccCutESD[200];
223   sprintf(ResonanceAccCutESD,"pt >=  %.2f && pt<= %.2f && minv>=  %.2f && minv<= %.2f",ptcutmin,ptcutmax,masscutmin,masscutmax);
224
225
226
227   /*********************************/
228   // Cut conditions (Id,Background...)
229
230   Char_t IdcutResonanceMC[100];      
231   Char_t IdcutMuminusMC[100];   
232   Char_t IdcutMuplusMC[100];      
233  
234   if (ResType==553){
235     sprintf(IdcutResonanceMC,"id==553");
236     sprintf(IdcutMuminusMC,"id==13 && idmo==553");
237     sprintf(IdcutMuplusMC,"id==-13 && idmo==553");  
238   }
239   if (ResType==443){
240     sprintf(IdcutResonanceMC,"id==443");
241     sprintf(IdcutMuminusMC,"id==13 && idmo==443");
242     sprintf(IdcutMuplusMC,"id==-13 && idmo==443");  
243   }
244   
245   //means no cuts since we don't have the trackid propagated yet pt>0
246   // now we have it 10/05/2005 but it's still not being used
247   
248   
249   // Background is calculated in the MUONefficiency.C macro 
250   // Background calculation is meaningful only when Resonance have been generated with realistic background
251   // when it's a pure Resonance file, the background lies artificially at the Resonance mass 
252   
253   //no realistic background
254   Int_t REALISTIC_BACKGROUND = 0 ; 
255
256   // to take background into account, i.e. substraction, if necessary
257   // this is not ready yet
258   Char_t BckgdCutResonanceESD[100];      
259   sprintf(BckgdCutResonanceESD,"pt>0 && minv>%f && minv<%f && pt1>0 && pt2>0",masscutmin,masscutmax); 
260   
261   // if realistic background 
262   // same cut + substract the background from hInvMassBg
263   
264
265
266   /*******************************************************************************************************************/
267   Char_t txt[50];
268   TLatex *tex;
269   TLegend *leg;
270   
271   /*******************************/
272   /*      Monte Carlo Part       */
273   /*******************************/
274
275   //----------------------------
276   // Pt-rapidity distributions from Kinematics
277   //----------------------------
278   TH2F *hptyResonanceMC = new TH2F("hptyResonanceMC", " Monte Carlo Resonance",ptbins,ptmin,ptmax,ybins,ymin,ymax);
279   hptyResonanceMC->SetLineColor(1);
280   hptyResonanceMC->SetLineStyle(1);
281   hptyResonanceMC->SetLineWidth(2);
282   Ktuple->Project("hptyResonanceMC","y:pt",IdcutResonanceMC);
283   hptyResonanceMC->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
284   hptyResonanceMC->GetYaxis()->SetTitle("Rapidity");
285  
286   cout << "  " << endl;
287   cout << "********************" << endl;
288   cout << " Monte Carlo Tracks  "<< endl;
289   cout << " " << hptyResonanceMC->GetEntries() << " Resonance in simulation " << endl;
290   
291
292   //******** Add acceptance cuts  - Theta and Pt
293   TH2F *hptyResonanceMCAcc = new TH2F("hptyResonanceMCAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax,ybins,ymin,ymax);
294   hptyResonanceMCAcc->SetLineColor(1);
295   hptyResonanceMCAcc->SetLineStyle(1);
296   hptyResonanceMCAcc->SetLineWidth(2);
297
298   TString m1MC(IdcutResonanceMC);
299   TString m2MC(ResonanceAccCutMC);
300   TString m3MC = m1MC + " && " + m2MC ;
301
302   Ktuple->Project("hptyResonanceMCAcc","y:pt",m3MC.Data());
303   hptyResonanceMCAcc->GetYaxis()->SetTitle("Rapidity");
304   hptyResonanceMCAcc->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
305   hptyResonanceMCAcc->SetStats(1);
306   hptyResonanceMCAcc->Sumw2();
307   
308   TCanvas *c1 = new TCanvas("c1", "Resonance  Monte Carlo:  Pt vs Y",30,30,700,500);
309   c1->Range(-1.69394,-0.648855,15.3928,2.77143);
310   c1->SetBorderSize(2);
311   c1->SetRightMargin(0.0229885);
312   c1->SetTopMargin(0.0275424);
313   c1->SetFrameFillColor(0);
314   c1->cd();
315
316   hptyResonanceMCAcc->SetStats(0);
317   hptyResonanceMCAcc->Draw("LEGO2ZFB"); 
318   //TLatex *tex = new TLatex(2,6,"#mu^{-}");
319   //tex->SetTextSize(0.06);
320   //tex->SetLineWidth(2);  
321   //tex->Draw();
322  
323   sprintf(txt,"Resonance : %d entries",(Int_t)hptyResonanceMCAcc->GetEntries());
324   tex = new TLatex(-0.854829,0.794436,txt);
325   tex->SetLineWidth(2);
326   tex->Draw();
327
328   c1->Update();
329   if (SAVE){
330     c1->SaveAs("ptyResonanceMCAcc.gif");
331     c1->SaveAs("ptyResonanceMCAcc.eps");
332   }
333
334
335   TH1F *hptResonanceMCAcc = new TH1F("hptResonanceMCAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax);
336   hptResonanceMCAcc->SetLineColor(1);
337   hptResonanceMCAcc->SetLineStyle(1);
338   hptResonanceMCAcc->SetLineWidth(2);
339   Ktuple->Project("hptResonanceMCAcc","pt",m3MC.Data());
340   hptResonanceMCAcc->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
341   hptResonanceMCAcc->SetStats(1);
342   hptResonanceMCAcc->Sumw2();
343   
344   TH1F *hyResonanceMCAcc = new TH1F("hyResonanceMCAcc", " Monte Carlo Resonance in acceptance",ybins,ymin,ymax);
345   hyResonanceMCAcc->SetLineColor(1);
346   hyResonanceMCAcc->SetLineStyle(1);
347   hyResonanceMCAcc->SetLineWidth(2);
348   Ktuple->Project("hyResonanceMCAcc","y",m3MC.Data());
349   hyResonanceMCAcc->GetXaxis()->SetTitle("Rapidity");
350   hyResonanceMCAcc->SetStats(1);
351   hyResonanceMCAcc->Sumw2();
352   
353   
354   //   TH2F *hKAccptyCutsMuplus = new TH2F("hKAccptyCutsMuplus", "Monte Carlo #mu^{+} ",ptbins,ptmin,ptmax,ybins,ymin,ymax);
355   //   hKAccptyCutsMuplus->SetLineColor(1);
356   //   hKAccptyCutsMuplus->SetLineStyle(1);
357   //   hKAccptyCutsMuplus->SetLineWidth(2);
358   
359   //   TString p1MC(IdcutMuplusMC);
360   //   TString p2MC(AccCutMC);
361   //   TString p3MC = p1MC + " && " + p2MC ;
362   
363   //   Ktuple->Project("hKAccptyCutsMuplus","y:pt",p3MC.Data());
364   //   hKAccptyCutsMuplus->GetYaxis()->SetTitle("Rapidity");
365   //   hKAccptyCutsMuplus->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
366   //   hKAccptyCutsMuplus->SetStats(1);
367   //   hKAccptyCutsMuplus->Sumw2();
368   
369   
370   cout << "  " << endl;
371   cout << "********************" << endl;
372   cout << " Monte Carlo Tracks  "<< endl;
373   cout << " " << hptyResonanceMCAcc->GetEntries() << " Resonance in acceptance cuts " << endl;
374   
375
376   /*******************************************************************************************************************/
377
378   /*******************************/
379   /*  Reconstructed Tracks Study */
380   /*******************************/
381
382   //----------------------------
383   // Pt-rapidity distributions from ESD : reconstructed tracks/particle
384   //----------------------------
385   TH2F *hptyResonanceESD = new TH2F("hptyResonanceESD", " Reconstucted Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax);
386   hptyResonanceESD->SetLineColor(1);
387   hptyResonanceESD->SetLineStyle(1);
388   hptyResonanceESD->SetLineWidth(2);
389   ESDtuple->Project("hptyResonanceESD","y:pt",BckgdCutResonanceESD);
390   hptyResonanceESD->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
391   hptyResonanceESD->GetYaxis()->SetTitle("Rapidity");
392
393   cout << "  " << endl;
394   cout << "********************" << endl;
395   cout << " Reconstructed Tracks" << endl ;
396   cout << " " << hptyResonanceESD->GetEntries() << " Resonance reconstructed " << endl;
397
398   if (REALISTIC_BACKGROUND){
399     TH2F *hptyResonanceESDBck = new TH2F("hptyResonanceESDBck", "Background",ptbins,ptmin,ptmax,ybins,ymin,ymax);
400     hptyResonanceESDBck->SetLineColor(1);
401     hptyResonanceESDBck->SetLineStyle(1);
402     hptyResonanceESDBck->SetLineWidth(2);
403     ESDtupleBck->Project("hptyResonanceESDBck","y:pt",BckgdCutResonanceESD);
404     hptyResonanceESDBck->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
405     hptyResonanceESDBck->GetYaxis()->SetTitle("Rapidity");
406     cout << " with " << hptyResonanceESDBck->GetEntries() << " Resonances from Background (random mixing) " << endl;
407   }
408
409   // if something is wrong 
410   if ( hptyResonanceESD->GetEntries()==0) {cout << " No entries in hptyResonanceESD " << endl ; return 0 ;}
411
412   //with Acc cuts - Theta and Pt
413   TH2F *hptyResonanceESDAcc = new TH2F("hptyResonanceESDAcc", "Reconstructed Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax);
414   hptyResonanceESDAcc->SetLineColor(1);
415   hptyResonanceESDAcc->SetLineStyle(1);
416   hptyResonanceESDAcc->SetLineWidth(2);
417   
418   TString m1Rec(BckgdCutResonanceESD);
419   TString m2Rec(ResonanceAccCutESD);
420   TString m3Rec = m1Rec + " && " + m2Rec ;
421
422   ESDtuple->Project("hptyResonanceESDAcc","y:pt",m3Rec.Data());
423   hptyResonanceESDAcc->GetYaxis()->SetTitle("Rapidity");
424   hptyResonanceESDAcc->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
425   hptyResonanceESDAcc->SetStats(1);
426   hptyResonanceESDAcc->Sumw2();
427
428   cout << "  " << endl;
429   cout << "********************" << endl;
430   cout << " Reconstructed Tracks" << endl ;
431   cout << " " << hptyResonanceESDAcc->GetEntries() << " Resonance in acceptance cuts " << endl;
432
433
434   TH2F *hptyResonanceESDBckAcc;
435   if (REALISTIC_BACKGROUND){
436     //with Acc cuts - Theta and Pt
437     hptyResonanceESDBckAcc = new TH2F("hptyResonanceESDBckAcc", "Reconstructed Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax);
438     hptyResonanceESDBckAcc->SetLineColor(1);
439     hptyResonanceESDBckAcc->SetLineStyle(1);
440     hptyResonanceESDBckAcc->SetLineWidth(2);
441     
442     TString m1Rec(BckgdCutResonanceESD);
443     TString m2Rec(ResonanceAccCutESD);
444     TString m3Rec = m1Rec + " && " + m2Rec ;
445     
446     ESDtupleBck->Project("hptyResonanceESDBckAcc","y:pt",m3Rec.Data());
447     hptyResonanceESDBckAcc->GetYaxis()->SetTitle("Rapidity");
448     hptyResonanceESDBckAcc->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
449     hptyResonanceESDBckAcc->SetStats(1);
450     hptyResonanceESDBckAcc->Sumw2();
451     cout << " with " << hptyResonanceESDBckAcc->GetEntries() << " Resonances from Background (random mixing) " << endl;
452   }
453
454
455   TCanvas *c100 = new TCanvas("c100", "Resonance Reconstructed in Acceptance: Pt vs Y",210,30,700,500);
456   c100->Range(-1.69394,-0.648855,15.3928,2.77143);
457   c100->SetBorderSize(2);
458   c100->SetRightMargin(0.0229885);
459   c100->SetTopMargin(0.0275424);
460   c100->SetFrameFillColor(0);
461   
462   c100->cd();
463   hptyResonanceESDAcc->SetStats(0);
464   hptyResonanceESDAcc->Draw("LEGO2ZFB"); 
465   sprintf(txt,"Resonance : %d entries",(Int_t)hptyResonanceESDAcc->GetEntries());
466   tex = new TLatex(-0.854829,0.794436,txt);
467   tex->SetLineWidth(2);
468   tex->Draw();
469
470   c100->Update();
471   if (SAVE){
472     c100->SaveAs("ptyResonanceESDAcc.gif");
473     c100->SaveAs("ptyResonanceESDAcc.eps");
474   }
475
476   if (REALISTIC_BACKGROUND){
477     TCanvas *c110 = new TCanvas("c110", "Resonance Background Reconstructed in Acceptance: Pt vs Y",215,35,700,500);
478     c110->Range(-1.69394,-0.648855,15.3928,2.77143);
479     c110->SetBorderSize(2);
480     c110->SetRightMargin(0.0229885);
481     c110->SetTopMargin(0.0275424);
482     c110->SetFrameFillColor(0);
483     
484     c110->cd();
485     hptyResonanceESDBckAcc->SetStats(0);
486     hptyResonanceESDBckAcc->Draw("LEGO2ZFB"); 
487     sprintf(txt,"Resonance backround : %d entries",(Int_t)hptyResonanceESDBckAcc->GetEntries());
488     tex = new TLatex(-0.854829,0.794436,txt);
489     tex->SetLineWidth(2);
490     tex->Draw();
491     
492     
493     c110->Update();
494     if (SAVE){
495       c110->SaveAs("ptyResonanceESDBckAcc.gif");
496       c110->SaveAs("ptyResonanceESDBckAcc.eps");
497     }
498   }
499   
500   TH1F *hptResonanceESDAcc = new TH1F("hptResonanceESDAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax);
501   hptResonanceESDAcc->SetLineColor(1);
502   hptResonanceESDAcc->SetLineStyle(1);
503   hptResonanceESDAcc->SetLineWidth(2);
504   ESDtuple->Project("hptResonanceESDAcc","pt",m3Rec.Data());
505   hptResonanceESDAcc->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
506   hptResonanceESDAcc->SetStats(1);
507   hptResonanceESDAcc->Sumw2();
508   
509   TH1F *hyResonanceESDAcc = new TH1F("hyResonanceESDAcc", " Monte Carlo Resonance in acceptance",ybins,ymin,ymax);
510   hyResonanceESDAcc->SetLineColor(1);
511   hyResonanceESDAcc->SetLineStyle(1);
512   hyResonanceESDAcc->SetLineWidth(2);
513   ESDtuple->Project("hyResonanceESDAcc","y",m3Rec.Data());
514   hyResonanceESDAcc->GetXaxis()->SetTitle("Rapidity");
515   hyResonanceESDAcc->SetStats(1);
516   hyResonanceESDAcc->Sumw2();
517   
518
519   /*******************************************************************************************************************/
520
521   /*******************************/
522   /* Efficiencies calculations   */
523   /*******************************/
524   
525
526   cout << "  " << endl;
527   cout << "***********************" << endl;
528   cout << " Integrated efficiency" << endl ;
529   
530   if (REALISTIC_BACKGROUND)
531     cout << " ResonanceESDAcc/ResonanceMCAcc = " << (hptyResonanceESDAcc->GetEntries()-hptyResonanceESDBckAcc->GetEntries())/hptyResonanceMCAcc->GetEntries()  << endl;
532   else 
533     cout << " ResonanceESDAcc/ResonanceMCAcc = " << hptyResonanceESDAcc->GetEntries()/hptyResonanceMCAcc->GetEntries()  << endl;
534
535   TH2F *hEffptyResonance = new TH2F("hEffptyResonance", " Resonance Efficiency",ptbins,ptmin,ptmax,ybins,ymin,ymax);
536
537   if (REALISTIC_BACKGROUND){
538     TH2F *hptyResonanceESDAccBckSubstracted = new TH2F("hptyResonanceESDAccBckSubstracted","hptyResonanceESDAccBckSubstracted",ptbins,ptmin,ptmax,ybins,ymin,ymax);
539     hptyResonanceESDAccBckSubstracted->Add(hptyResonanceESDAcc,hptyResonanceESDBckAcc,1,-1);
540     hEffptyResonance->Divide(hptyResonanceESDAccBckSubstracted,hptyResonanceMCAcc,1,1);
541   }
542   else 
543   hEffptyResonance->Divide(hptyResonanceESDAcc,hptyResonanceMCAcc,1,1);
544
545   TCanvas *c1000 = new TCanvas("c1000", "Resonance efficiency : Pt vs Y",390,30,700,500);
546   c1000->Range(-1.69394,-0.648855,15.3928,2.77143);
547   c1000->SetBorderSize(2);
548   c1000->SetRightMargin(0.0229885);
549   c1000->SetTopMargin(0.0275424);
550   c1000->SetFrameFillColor(0);
551
552   c1000->cd();
553   hEffptyResonance->SetStats(0);
554   hEffptyResonance->Draw("LEGO2fz");
555
556
557
558   TH1F *hEffptResonance = new TH1F("hEffptResonance", "Resonance Efficiency vs pt",ptbins,ptmin,ptmax);
559   hEffptResonance->Divide(hptResonanceESDAcc,hptResonanceMCAcc,1,1);
560   hEffptResonance->SetLineWidth(2);
561   hEffptResonance->SetMinimum(0);
562   hEffptResonance->SetStats(1);
563
564
565   TCanvas *c1100 = new TCanvas("c1100", "Resonance efficiency : Pt",410,50,700,500);
566   c1100->Range(-1.69394,-0.648855,15.3928,2.77143);
567   c1100->SetBorderSize(2);
568   c1100->SetRightMargin(0.0229885);
569   c1100->SetTopMargin(0.0275424);
570   c1100->SetFrameFillColor(0);
571
572   c1100->cd();
573   hEffptResonance->SetStats(0);
574   hEffptResonance->GetXaxis()->SetTitle("P_{#perp}  [GeV/c]");
575   hEffptResonance->GetYaxis()->SetTitle("Efficiency ");
576   hEffptResonance->Draw("E");
577    
578   
579   
580   TH1F *hEffyResonance = new TH1F("hEffyResonance", "Resonance Efficiency vs y",ybins,ymin,ymax);
581   hEffyResonance->Divide(hyResonanceESDAcc,hyResonanceMCAcc,1,1);
582   hEffyResonance->SetLineWidth(2);
583   hEffyResonance->SetStats(1);
584   
585   TCanvas *c1200 = new TCanvas("c1200", "Resonance efficiency : Y",430,70,700,500);
586   c1200->Range(-1.69394,-0.648855,15.3928,2.77143);
587   c1200->SetBorderSize(2);
588   c1200->SetRightMargin(0.0229885);
589   c1200->SetTopMargin(0.0275424);
590   c1200->SetFrameFillColor(0);
591
592   c1200->cd();
593   hEffyResonance->SetStats(0);
594   hEffyResonance->GetXaxis()->SetTitle("Rapidity");
595   hEffyResonance->GetYaxis()->SetTitle("Efficiency ");
596   hEffyResonance->Draw("E");
597
598   c1000->Update();
599   c1100->Update();
600   c1200->Update();
601   if (SAVE){
602     c1000->SaveAs("EffptyResonance.gif");
603     c1000->SaveAs("EffptyResonance.eps");
604     c1100->SaveAs("EffptResonance.gif");
605     c1100->SaveAs("EffptResonance.eps");
606     c1200->SaveAs("EffyResonance.gif");
607     c1200->SaveAs("EffyResonance.eps");
608   }
609
610  /*******************************************************************************************************************/
611
612   /*******************************************/
613   /* Trigger matching                        */
614   /* trigger word is defined in :            */
615   /* /STEER/createTriggerDescriptor_MUON.C   */
616   /*******************************************/
617
618
619   // Float_t triggerChi2Min = 0.; 
620   // Float_t triggerChi2Max = 7.5;
621
622   //TString m1Rec(BckgdCutResonanceESD);
623   //TString m2Rec(ResonanceAccCutESD);
624   TString  m1TG = m1Rec + " && " + m2Rec ;
625   //cout << m1TG.Data() << endl;
626   
627   TH1F *hInvMassNoTriggerCut= new TH1F("hInvMassNoTriggerCut","hInvMassNoTriggerCut",invMassBins,invMassMin,invMassMax);
628   hInvMassNoTriggerCut->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
629   hInvMassNoTriggerCut->GetYaxis()->SetTitle("Counts");
630   ESDtuple->Project("hInvMassNoTriggerCut","minv",m1TG.Data());
631
632   //Add trigger UnlikeHighPt or UnlikeLowPt
633   TString  m2TG = m1TG + " && ((tw & 0x10) == 16 || (tw & 0x32) == 32)" ;
634   TH1F *hInvMassResonanceTrigger= new TH1F("hInvMassResonanceTrigger","hInvMassResonanceTrigger",invMassBins,invMassMin,invMassMax);
635   hInvMassResonanceTrigger->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
636   hInvMassResonanceTrigger->GetYaxis()->SetTitle("Counts");
637   ESDtuple->Project("hInvMassResonanceTrigger","minv",m2TG.Data());
638
639
640   cout << "  " << endl;
641   cout << "********************" << endl;
642   cout << " TriggerMatching" << endl ;
643   cout << " " << hInvMassResonanceTrigger->GetEntries() << " Resonance with trigger UnlikePairAllPt " << endl;
644
645   TString  m2TGNo = m1TG + " && (tw & 0x800) != 2048" ;
646   TH1F *hResonanceTriggerNo= new TH1F("hResonanceTriggerNo","hResonanceTriggerNo",32769,0,32768);
647   hResonanceTriggerNo->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
648   hResonanceTriggerNo->GetYaxis()->SetTitle("Counts");
649   ESDtuple->Project("hResonanceTriggerNo","tw",m2TGNo.Data());
650
651   
652
653
654   //Add matching rec/trig for 2 tracks 
655   TString m3TG = m2TG + " && trig1 > 0 && trig2 > 0" ;  
656   TH1F *hInvMassResonanceTriggerTwoMatch = new TH1F("hInvMassResonanceTriggerTwoMatch","hInvMassResonanceTrigger with 2 matching tracks",invMassBins,invMassMin,invMassMax);
657   hInvMassResonanceTriggerTwoMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
658   hInvMassResonanceTriggerTwoMatch->GetYaxis()->SetTitle("Counts");
659   ESDtuple->Project("hInvMassResonanceTriggerTwoMatch","minv",m3TG.Data());
660
661
662   //Add matching rec/trig for 1  tracks 
663   TString m4TG = m2TG + " && (trig1 > 0 || trig2 > 0)" ;  
664   TH1F *hInvMassResonanceTriggerOneMatch= new TH1F("hInvMassResonanceTriggerOneMatch","hInvMassResonanceTrigger with 1 matching tracks",invMassBins,invMassMin,invMassMax);
665   hInvMassResonanceTriggerOneMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
666   hInvMassResonanceTriggerOneMatch->GetYaxis()->SetTitle("Counts");
667   ESDtuple->Project("hInvMassResonanceTriggerOneMatch","minv",m4TG.Data());
668
669   TString m5TG = m2TG + " && (trig1 == 0 && trig2 == 0)" ;  
670   TH1F *hInvMassResonanceTriggerNoMatch= new TH1F("hInvMassResonanceTriggerNoMatch","hInvMassResonanceTrigger with no matching tracks",invMassBins,invMassMin,invMassMax);
671   hInvMassResonanceTriggerNoMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
672   hInvMassResonanceTriggerNoMatch->GetYaxis()->SetTitle("Counts");
673   ESDtuple->Project("hInvMassResonanceTriggerNoMatch","minv",m5TG.Data());
674
675   TCanvas *c2 = new TCanvas("c2", "Resonance Trigger efficiency",30,90,700,500);
676   c2->Range(-1.69394,-0.648855,15.3928,2.77143);
677   c2->SetBorderSize(2);
678   c2->SetRightMargin(0.0229885);
679   c2->SetTopMargin(0.0275424);
680   c2->SetFrameFillColor(0);
681   c2->SetLogy(1);
682
683   c2->cd();
684   hInvMassNoTriggerCut->SetStats(0);
685   hInvMassNoTriggerCut->SetLineWidth(2);
686   hInvMassNoTriggerCut->SetLineStyle(3);
687   hInvMassNoTriggerCut->SetLineColor(1);
688   hInvMassNoTriggerCut->GetYaxis()->SetTitleOffset(1.2);
689   hInvMassNoTriggerCut->Draw();
690   hInvMassResonanceTrigger->SetLineWidth(2);
691   hInvMassResonanceTrigger->SetLineStyle(0);
692   hInvMassResonanceTrigger->SetLineColor(2);
693   hInvMassResonanceTrigger->Draw("same");
694   hInvMassResonanceTriggerTwoMatch->SetLineWidth(2);
695   hInvMassResonanceTriggerTwoMatch->SetLineStyle(1);
696   hInvMassResonanceTriggerTwoMatch->SetLineColor(4);
697   hInvMassResonanceTriggerTwoMatch->Draw("same");
698   hInvMassResonanceTriggerOneMatch->SetLineWidth(2);
699   hInvMassResonanceTriggerOneMatch->SetLineStyle(2);
700   hInvMassResonanceTriggerOneMatch->SetLineColor(51);
701   hInvMassResonanceTriggerOneMatch->Draw("same");
702   hInvMassResonanceTriggerNoMatch->SetLineWidth(2);
703   hInvMassResonanceTriggerNoMatch->SetLineStyle(2);
704   hInvMassResonanceTriggerNoMatch->SetLineColor(46);
705   hInvMassResonanceTriggerNoMatch->Draw("same");
706
707   leg = new TLegend(0.12,0.6,0.50,0.89);
708   leg->SetHeader("Reconstructed Resonance Invariant Mass");
709
710   leg->AddEntry(hInvMassNoTriggerCut,Form("All  (%.0f cnts)",hInvMassNoTriggerCut->GetEntries()),"l");
711   leg->AddEntry(hInvMassResonanceTrigger,Form("UnlikePairAllPt Trigger (%.0f cnts)",hInvMassResonanceTrigger->GetEntries()),"l");
712   leg->AddEntry(hInvMassResonanceTriggerTwoMatch,Form("UPAllPt Trig. and 2 tracks matches (%.0f cnts)",hInvMassResonanceTriggerTwoMatch->GetEntries()),"l");  
713   leg->AddEntry(hInvMassResonanceTriggerOneMatch,Form("UPAllPt Trig. and 1 track match (%.0f cnts)",hInvMassResonanceTriggerOneMatch->GetEntries()),"l");  
714   leg->AddEntry(hInvMassResonanceTriggerNoMatch,Form("UPAllPt Trig. and no matching track (%.0f cnts)",hInvMassResonanceTriggerNoMatch->GetEntries()),"l");  
715   leg->Draw();
716
717   c2->Update();
718   if(SAVE){
719     c2->SaveAs("TriggerMatching.gif");
720     c2->SaveAs("TriggerMatching.eps");
721   }
722   
723  /*******************************************************************************************************************/
724
725   /*******************************/
726   /* Mass resolution             */
727   /*******************************/
728   
729   const Int_t nofMassHistograms = ptbins ; 
730   // cs invMassMin = 7.0 ;
731   // cs invMassMax = 11.0 ;
732   // cs invMassBins = 100 ;
733   
734   //Float_t ptmin = 0.0;      
735   //Float_t ptmax =  16.0;  Int_t ptbins = 8;    
736   Float_t ptbinssize = ptmax/ptbins; 
737   
738
739   Float_t ptbinslimits[nofMassHistograms+1];
740   Float_t ptbinscenters[nofMassHistograms];
741   
742   for(Int_t as = 0 ; as < nofMassHistograms+1 ; as++){
743     ptbinslimits[as] = as*ptbinssize ;
744     //cout << "ptbinslimits["<<as<<"] = " << ptbinslimits[as] << endl ;
745   }
746   for(Int_t as = 0 ; as < nofMassHistograms ; as++) {
747     ptbinscenters[as] = ptbinslimits[as] + (ptbinslimits[as+1] - ptbinslimits[as])/2 ;
748     //cout << "ptbinscenters[as] = " << ptbinscenters[as] << endl ;
749   }
750   
751
752   TH1F **hInvMassInPtBins = new TH1F* [nofMassHistograms] ;
753   TString histExt;
754   Char_t  theExt[20];
755   TString histCut;
756   Char_t  theCut[100];
757   
758   
759   TF1 *fitFunc = 0;  
760   TString FitFuncName;
761
762   if (fitfunc==0 || fitfunc==1){
763     fitFunc = new TF1("gaussian","gaus(0)",FitLow,FitHigh);
764     if (!fitFunc)
765       fitFunc = new TF1("gaussian","[0]*TMath::Exp(-0.5*(x-[1])*(x-[1])/[2]/[2])",FitLow,FitHigh);
766     fitFunc->SetParNames("Constant","Mean value","Width");
767     fitFunc->SetLineColor(2);
768     fitFunc->SetLineWidth(2);
769     fitFunc->SetLineStyle(2);
770     
771     Float_t  *tParams = new Float_t[3] ;
772     tParams[0] = 500;
773     tParams[1] = RESONANCE_MASS;
774     tParams[2] = 0.1;
775
776     fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]);
777     
778     FitFuncName = "gaussian";  
779   }
780   
781   if (fitfunc==2){
782     fitFunc = new TF1("fitlandau_a",fitlandau_a,FitLow,FitHigh,3);
783     fitFunc->SetParNames("Constant","Mean value","Width");
784     fitFunc->SetLineColor(2);
785     fitFunc->SetLineWidth(2);
786     fitFunc->SetLineStyle(2);
787     
788     Float_t  *tParams = new Float_t[3] ;
789     tParams[0] = 1000;
790     tParams[1] = RESONANCE_MASS;
791     tParams[2] = 0.05;
792     fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]);
793     
794     TString FitFuncName = "fitlandau_a";  
795   }
796   
797   if (fitfunc==3){
798     fitFunc = new TF1("fitlandaugauss",fitlandaugauss,FitLow,FitHigh,6);
799     fitFunc->SetParNames("Constant","Mean value","Width","SigmaGauss","LowFitVal","HighFitVal");
800     fitFunc->SetLineColor(2);
801     fitFunc->SetLineWidth(2);
802     fitFunc->SetLineStyle(2);
803     
804     Float_t  *tParams = new Float_t[6] ;
805     tParams[0] = 5000;
806     tParams[1] = RESONANCE_MASS;   
807     tParams[2] = 0.05;   //0.5
808     tParams[3] = 0.05;   // 1.
809     
810     tParams[4] = FitLow;
811     tParams[5] = FitHigh;
812
813     fitFunc->SetParameters(tParams[0],tParams[1],tParams[2],tParams[3],tParams[4],tParams[5]);
814     fitFunc->FixParameter(4,FitLow);
815     fitFunc->FixParameter(5,FitHigh);
816     
817     fitFunc->SetParLimits(1, 9.1 ,9.7);
818     fitFunc->SetParLimits(2, 0.001 ,0.1);
819     fitFunc->SetParLimits(3, 0.001 ,0.5);
820
821     // for special case one may use 
822     //fitFunc->SetParameter(1,9.35);
823
824     TString FitFuncName = "fitlandaugauss"; 
825   }
826   
827   if (fitfunc==4){
828     fitFunc = new TF1("fitlandau",fitlandau,FitLow,FitHigh,3);
829     fitFunc->SetParNames("Constant","Mean value","Width");
830     fitFunc->SetLineColor(2);
831     fitFunc->SetLineWidth(2);
832     fitFunc->SetLineStyle(2);
833     
834     Float_t  *tParams = new Float_t[3] ;
835     tParams[0] = 1000;
836     tParams[1] = RESONANCE_MASS;
837     tParams[2] = 0.05;
838     fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]);
839     
840     TString FitFuncName = "fitlandau";  
841   }
842
843
844   Float_t fitResultsWidth[nofMassHistograms];
845   Float_t fitResultsWidthErr[nofMassHistograms];
846   
847   Float_t fitResultsMean[nofMassHistograms];
848   Float_t fitResultsMeanErr[nofMassHistograms];
849   
850   Float_t errx[nofMassHistograms];
851
852
853
854   TCanvas *call = new TCanvas("call", "Resonance : invariant mass spectra",30,330,700,500);
855   call->Range(-1.69394,-0.648855,15.3928,2.77143);
856   call->SetBorderSize(2);
857   //call->SetRightMargin(0.2229885);
858   call->SetRightMargin(0.0229885);
859   call->SetTopMargin(0.0275424);
860   call->SetFrameFillColor(0);
861   call->cd();
862
863   TH1F* hInvMass = new TH1F("hInvMass","Inv. Mass Pt,Y integrated",30,0,3+RESONANCE_MASS);
864   ESDtuple->Project("hInvMass","minv");
865   hInvMass->SetLineWidth(2);
866   hInvMass->Draw("HE");
867
868   if(REALISTIC_BACKGROUND){  
869     TH1F* hInvMassBck = new TH1F("hInvMassBck","Background Pt,Y integrated",30,0,3+RESONANCE_MASS);
870     ESDtupleBck->Project("hInvMassBck","minv");
871     hInvMassBck->SetLineWidth(2);
872     hInvMassBck->SetLineStyle(2);
873     hInvMassBck->SetLineColor(2);    
874     hInvMassBck->Draw("same");
875   }
876   call->Modified();
877   call->Update();
878
879
880   TCanvas *cfit = new TCanvas("cfit", "Resonance : invariant mass fit",30,330,700,500);
881   cfit->Range(-1.69394,-0.648855,15.3928,2.77143);
882   cfit->SetBorderSize(2);
883   //cfit->SetRightMargin(0.2229885);
884   cfit->SetRightMargin(0.0229885);
885   cfit->SetTopMargin(0.0275424);
886   cfit->SetFrameFillColor(0);
887   
888   //TH1F *hInvMassAllClone = (TH1F*) gROOT->FindObject("hInvMassAll");
889   
890   TH1F* hInvMassAll = new TH1F("hInvMassAll","Inv. Mass Pt,Y integrated",invMassBins,invMassMin,invMassMax);
891   sprintf(theCut,"pt> %.2f && pt<= %.2f",ptbinslimits[0],ptbinslimits[nofMassHistograms]);
892   //cout << "theCut" << theCut << endl ;
893   ESDtuple->Project("hInvMassAll","minv",theCut);
894   hInvMassAll->Draw();
895   
896   cfit->Update();
897   if (SAVE){
898     cfit->SaveAs("ResonanceMass.gif");
899     cfit->SaveAs("ResonanceMass.eps");
900   }
901
902   if(REALISTIC_BACKGROUND){  
903     TH1F* hInvMassAllBck = new TH1F("hInvMassAllBck","Background Pt,Y integrated",invMassBins,invMassMin,invMassMax);
904     ESDtupleBck->Project("hInvMassAllBck","minv",theCut);
905     hInvMassAllBck->SetLineWidth(2);
906     hInvMassAllBck->SetLineStyle(2);
907     hInvMassAllBck->SetLineColor(2);    
908     hInvMassAllBck->Draw("same");
909   }
910   cfit->Modified();
911   cfit->Update();
912   
913   //Fit also the pt-integrated mass histogram 
914   cfit->cd();
915   hInvMassAll->Fit(FitFuncName.Data(),"rv");
916   if (FitFuncName == "gaussian")
917     fitFunc->SetRange(fitFunc->GetParameter(1)-gaussianFitWidth,fitFunc->GetParameter(1)+gaussianFitWidth);
918   hInvMassAll->Fit(FitFuncName.Data(),"rv");
919   cfit->Modified();
920   cfit->Update();
921   
922   cout << " " << endl;
923   cout << "********************************************" << endl;
924   cout << " Fit ("<<FitFuncName.Data()<<" results of the pt-integrated mass peak " << endl ;
925    cout << "   Mean value : " << fitFunc->GetParameter(1) << " +/- " << fitFunc->GetParError(1) << endl ;
926    cout << "   Width value : " << fitFunc->GetParameter(2) << " +/- " << fitFunc->GetParError(2) << endl ;
927    cout << "   ChiSquare of the fit : " << fitFunc->GetChisquare()<< " / " << fitFunc->GetNDF() << endl ;
928    cout << "********************************************" << endl;
929    cout << " " << endl;
930    
931    TCanvas *cptfits = new TCanvas("cptfits", "Resonance : invariant mass fit",30,330,700,500);
932    cptfits->Range(-1.69394,-0.648855,15.3928,2.77143);
933    cptfits->SetBorderSize(2);
934    //cptfits->SetRightMargin(0.2229885);
935    cptfits->SetRightMargin(0.0229885);
936    cptfits->SetTopMargin(0.0275424);
937    cptfits->SetFrameFillColor(0);
938
939    Float_t chi2sum = 0.;
940   
941   for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++){
942     sprintf(theExt,"_%1f_%1f",ptbinslimits[qw],ptbinslimits[qw+1]);
943     histExt= theExt;
944     hInvMassInPtBins[qw] = new TH1F("hInvMassInPtBins"+histExt,"hInvMassInPtBins"+histExt,invMassBins,invMassMin,invMassMax);
945     hInvMassInPtBins[qw]->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");
946     hInvMassInPtBins[qw]->GetYaxis()->SetTitle("Counts");
947      
948     sprintf(theCut,"pt> %.2f && pt<= %.2f",ptbinslimits[qw],ptbinslimits[qw+1]);
949     histCut = theCut; 
950     ESDtuple->Project("hInvMassInPtBins"+histExt,"minv",histCut.Data());
951     
952     if (hInvMassInPtBins[qw]->GetEntries() > 80){
953       cptfits->cd();
954
955       if (FitFuncName == "gaussian"){
956         //prefit with a gaussian to get the maximum set a correct range
957         hInvMassInPtBins[qw]->Fit(FitFuncName.Data(),"rv");
958         fitFunc->SetRange(fitFunc->GetParameter(1)-gaussianFitWidth,fitFunc->GetParameter(1)+gaussianFitWidth);
959       }
960       
961       hInvMassInPtBins[qw]->SetStats(1);
962       hInvMassInPtBins[qw]->Draw("HE");
963       hInvMassInPtBins[qw]->Fit(FitFuncName.Data(),"rv");
964
965       
966       cout << "==============" << endl;
967       cout << "ChiSquare of the fit : " << fitFunc->GetChisquare()<< " / " << fitFunc->GetNDF() 
968            << " = " << fitFunc->GetChisquare()/fitFunc->GetNDF() <<endl;
969       cout << "==============" << endl;
970       chi2sum += fitFunc->GetChisquare()/fitFunc->GetNDF();
971  
972       cptfits->Modified();
973       cptfits->Update();
974       
975       fitResultsWidth[qw] = TMath::Abs(fitFunc->GetParameter(2));
976       fitResultsWidthErr[qw] =  fitFunc->GetParError(2);
977
978       fitResultsMean[qw] = TMath::Abs(fitFunc->GetParameter(1));
979       fitResultsMeanErr[qw] =  fitFunc->GetParError(1);
980       
981       if(FitFuncName == "fitlandaugauss"){      
982         // width = gauss_width + width landau
983         fitResultsWidth[qw] = TMath::Abs(fitFunc->GetParameter(2)) +  TMath::Abs(fitFunc->GetParameter(3)); 
984         fitResultsWidthErr[qw] =  TMath::Abs(fitFunc->GetParError(2)) +  TMath::Abs(fitFunc->GetParError(3));
985         
986         // problem with the mean of the fit (parameter1)  which doesn't give the peak position 
987         fitResultsMean[qw] = TMath::Abs(fitFunc->GetMaximumX());
988         fitResultsMeanErr[qw] =  fitFunc->GetParError(1); 
989         //this could work but it's not very nice....
990       }
991
992     }
993     else {
994       fitResultsWidth[qw] = 0.;
995       fitResultsWidthErr[qw] = 0.;
996       fitResultsMean[qw] = 0.;
997       fitResultsMeanErr[qw] = 0.;
998     }
999
1000     errx[qw] = 0.;
1001     
1002   }
1003   
1004   
1005   //From fit results extract a Mean value for the Mass and width
1006   cout << " " << endl;
1007   cout << "********************************************" << endl;
1008   cout << " Fit Function : " << FitFuncName  << " " ; 
1009   if (FitFuncName == "gaussian" ) cout << "+/- " << gaussianFitWidth << " MeV/c2" ;  
1010   cout <<  endl ;  
1011   cout << " Fit results :  mean values  " << endl ;
1012   Float_t meanMass = 0 ;
1013   Float_t meanMassError = 0 ;
1014   Float_t meanWidth = 0 ;
1015   Float_t meanWidthError = 0 ;
1016   Int_t cnt = 0 ;
1017   for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++){
1018     if ( fitResultsMean[qw] > invMassMin &&   fitResultsMean[qw] < invMassMax &&  fitResultsWidth[qw] > 0 &&  fitResultsWidth[qw] < 1){
1019       meanWidth += fitResultsWidth[qw] ;
1020       meanWidthError += fitResultsWidthErr[qw];
1021       meanMass += fitResultsMean[qw] ;
1022       meanMassError += fitResultsMeanErr[qw];
1023       cnt++;
1024     }
1025   }
1026   if (cnt==0) {
1027     cout << "Fitting procedure didn't work" << endl;  
1028     return 0 ;
1029   }
1030
1031   cout << " Mass : " << meanMass/cnt << " +/- " << meanMassError/cnt << endl ;
1032   cout << " Width : " << meanWidth/cnt << " +/- " << meanWidthError/cnt << endl ;
1033   cout << " Mean Chi2 of the fits : " << chi2sum/nofMassHistograms << endl ;
1034   cout << "********************************************" << endl;
1035   cout << " " << endl;
1036   
1037
1038   cout << "  " << endl;
1039   cout << "***********************" << endl;
1040   cout << " Integrated efficiency (+/- "<< 3*masssigma << " GeV around Resonance mass cut)" << endl ;
1041   cout << " ResonanceESDAcc/ResonanceMCAcc = " << hptyResonanceESDAcc->GetEntries() << "/" << hptyResonanceMCAcc->GetEntries()  << " = "  << hptyResonanceESDAcc->GetEntries()/hptyResonanceMCAcc->GetEntries()  <<endl;
1042   cout << "***********************" << endl;
1043   cout << "  " << endl;
1044
1045
1046   cout << "  " << endl;
1047   cout << "***********************" << endl;
1048   cout << " Trigger  efficiency" << endl ;
1049   cout << " Two muons matching = " << hInvMassResonanceTriggerTwoMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " <<  hInvMassResonanceTriggerTwoMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() << endl;
1050   cout << " Single muon matching  = " << hInvMassResonanceTriggerOneMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " <<  hInvMassResonanceTriggerOneMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() << endl;
1051   cout << " No matching  = " << hInvMassResonanceTriggerNoMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " <<  hInvMassResonanceTriggerNoMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() <<endl;
1052   cout << "***********************" << endl;
1053   cout << "  " << endl;
1054
1055   
1056   
1057   TCanvas *cA = new TCanvas("cA", "Resonance : Width from fit",30,330,700,500);
1058   cA->Range(-1.69394,-0.648855,15.3928,2.77143);
1059   cA->SetBorderSize(2);
1060   cA->SetRightMargin(0.0229885);
1061   cA->SetTopMargin(0.0275424);
1062   cA->SetFrameFillColor(0);
1063   cA->cd();
1064   
1065   //TH2F *hFrameA = new TH2F("hFrameA", "A hFrameA",ptbins,ptmin,ptmax,100,fitResultsWidth[1]-fitResultsWidth[1]*2,fitResultsWidth[1]+fitResultsWidth[1]*2);
1066   TH2F *hFrameA = new TH2F("hFrameA", "A hFrameA",ptbins,ptmin,ptmax,100,0.,2*fitResultsWidth[1]);
1067   hFrameA->GetYaxis()->SetTitle("Peak Width [GeV/c^{2}]");
1068   hFrameA->GetYaxis()->SetTitleOffset(1.2);
1069   hFrameA->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
1070   hFrameA->SetStats(0);
1071   hFrameA->Draw(); 
1072   
1073   
1074   TGraphErrors *grWidth = new TGraphErrors(nofMassHistograms,ptbinscenters,fitResultsWidth,errx,fitResultsWidthErr);
1075   grWidth->SetTitle("Resonance Mass Width");
1076   grWidth->SetMarkerColor(4);
1077   grWidth->SetMarkerStyle(21);
1078   grWidth->Draw("P");
1079   
1080   
1081   
1082   TCanvas *cB = new TCanvas("cB", "Resonance : Mean from fit",50,350,700,500);
1083   cB->Range(-1.69394,-0.648855,15.3928,2.77143);
1084   cB->SetBorderSize(2);
1085   cB->SetRightMargin(0.0229885);
1086   cB->SetTopMargin(0.0275424);
1087   cB->SetFrameFillColor(0);
1088   cB->cd();
1089   
1090   //TH2F *hFrameB = new TH2F("hFrameB", "A hFrameB",ptbins,ptmin,ptmax,100,fitResultsMean[1]-fitResultsMean[1]*2,fitResultsMean[1]+fitResultsMean[1]*2);
1091   TH2F *hFrameB = new TH2F("hFrameB", "A hFrameB",ptbins,ptmin,ptmax,100,invMassMin,invMassMax);
1092   hFrameB->GetYaxis()->SetTitle("Peak Position [GeV/c^{2}]");
1093   hFrameB->GetXaxis()->SetTitle("P_{#perp}   [GeV/c]");
1094   hFrameB->SetStats(0);
1095   hFrameB->Draw(); 
1096   
1097   
1098   
1099   
1100   //  TGaxis *Mean = new TGaxis(16,fitResultsWidth[1]-fitResultsWidth[1]*2,16,fitResultsWidth[1]+fitResultsWidth[1]*2,7,12,510,"+");
1101   //  Mean->SetLabelOffset(0.050);
1102   //   Mean->SetLabelSize(0.04);
1103   //   Mean->SetTickSize(0.03);
1104   //   Mean->SetGridLength(0);
1105   //   Mean->SetTitleOffset(0.1);
1106   //   Mean->SetTitleOffset(-1.1);
1107   //   Mean->SetTitleSize(0.05);
1108   //   Mean->SetTitle("Mean [Gev/c^2]");
1109   //   Mean->Draw();
1110   
1111   
1112   TGraphErrors *grMean = new TGraphErrors(nofMassHistograms,ptbinscenters,fitResultsMean,errx,fitResultsMeanErr);
1113   grMean->SetTitle("Resonance Mass Mean");
1114   grMean->SetMarkerColor(4);
1115   grMean->SetMarkerStyle(21);
1116   grMean->Draw("P");
1117   
1118   cA->Update();
1119   cB->Update();
1120   if (SAVE){
1121     cB->SaveAs("ResonanceMeanVsPt.gif");
1122     cB->SaveAs("ResonanceMeanVsPt.eps");
1123     cA->SaveAs("ResonanceWidthVsPt.gif");
1124     cA->SaveAs("ResonanceWidthVsPt.eps");
1125   }
1126   
1127   if (WRITE){
1128     TFile *myFile = new TFile("MUONplotefficiency.root", "RECREATE");
1129     hptyResonanceMCAcc->Write();
1130     hptyResonanceMC->Write();
1131     hptResonanceMCAcc ->Write();
1132     hyResonanceMCAcc->Write();
1133
1134     hptyResonanceESDAcc->Write();
1135     hptyResonanceESD->Write();
1136     hptResonanceESDAcc ->Write();
1137     hyResonanceESDAcc->Write();
1138
1139     hEffptyResonance->Write();
1140     hEffyResonance->Write();
1141     hEffptResonance->Write();
1142
1143     hInvMass->Write();
1144     hInvMassResonanceTrigger->Write();
1145     hResonanceTriggerNo->Write();
1146     hInvMassResonanceTriggerOneMatch->Write();
1147     hInvMassResonanceTriggerTwoMatch->Write();
1148     hInvMassResonanceTriggerNoMatch->Write();
1149
1150     for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++) 
1151       hInvMassInPtBins[qw]->Write();
1152     
1153     hInvMassAll->Write();
1154     
1155     myFile->Close();
1156   }
1157
1158   return 1;
1159   
1160 } // macro ends
1161
1162