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