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