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