PWG/muon/AliAnalysisTaskMuonDistributions.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 //-----------------------------------------------------------------------------
  19 // Analysis task to compute muon/dimuon kinematic distributions
  20 // The output is a list of histograms.
  21 // The macro class can run on AOD or in the train with the ESD filter.
  22 // R. Arnaldi
  23 //
  24 //-----------------------------------------------------------------------------
  25
  26 //#ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
  27 //#define ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
  28
  29 #include <TCanvas.h>
  30 #include <TF1.h>
  31 #include <TH1.h>
  32 #include <TLatex.h>
  33 #include <TList.h>
  34 #include <TStyle.h>
  35
  36 #include "AliAnalysisTaskMuonDistributions.h"
  37 #include "AliAnalysisTaskSE.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESD.h"
  41 #include "AliVEvent.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliInputEventHandler.h"
  44 #include "AliMCEvent.h"
  45 #include "AliStack.h"
  46 #include "AliLog.h"
  47 #include "AliHeader.h"
  48 #include "AliESDHeader.h"
  49 #include "AliStack.h"
  50 #include "TParticle.h"
  51 #include "TLorentzVector.h"
  52 #include "AliESDMuonTrack.h"
  53 #include "AliESDInputHandler.h"
  54 #include "AliAODEvent.h"
  55 #include "AliAODHeader.h"
  56 #include "AliAODHandler.h"
  57 #include "AliAODInputHandler.h"
  58
  59 ClassImp(AliAnalysisTaskMuonDistributions)
  60
  61 //__________________________________________________________________________
  62 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions() :
  63   fBeamEnergy(0.),
  64   fInvMassFitLimitMin(2.),
  65   fInvMassFitLimitMax(5.),
  66   fPsiFitLimitMin(2.9),
  67   fPsiFitLimitMax(3.3),
  68   fPsiPFitLimitMin(3.3),
  69   fPsiPFitLimitMax(4.),
  70   fBckFitLimitMin(2.2),
  71   fBckFitLimitMax(2.85),
  72   fInvariantMassFit(kFALSE),
  73   fkAnalysisType(0x0),
  74   fOutput(0x0)
  75 {
  76 }
  77 //___________________________________________________________________________
  78 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const Char_t* name) :
  79   AliAnalysisTaskSE(name),
  80   fBeamEnergy(0.),
  81   fInvMassFitLimitMin(2.),
  82   fInvMassFitLimitMax(5.),
  83   fPsiFitLimitMin(2.9),
  84   fPsiFitLimitMax(3.3),
  85   fPsiPFitLimitMin(3.3),
  86   fPsiPFitLimitMax(4.),
  87   fBckFitLimitMin(2.2),
  88   fBckFitLimitMax(2.85),
  89   fInvariantMassFit(kFALSE),
  90   fkAnalysisType(0x0),
  91   fOutput(0x0)
  92 {
  93   //
  94   // Constructor. Initialization of Inputs and Outputs
  95   //
  96   Info("AliAnalysisTaskMuonDistributions","Calling Constructor");
  97
  98   DefineOutput(1,TList::Class());
  99
 100 }
 101
 102 //___________________________________________________________________________
 103 AliAnalysisTaskMuonDistributions& AliAnalysisTaskMuonDistributions::operator=(const AliAnalysisTaskMuonDistributions& c)
 104 {
 105   //
 106   // Assignment operator
 107   //
 108   if (this!=&c) {
 109     AliAnalysisTaskSE::operator=(c) ;
 110   }
 111   return *this;
 112 }
 113
 114 //___________________________________________________________________________
 115 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const AliAnalysisTaskMuonDistributions& c) :
 116   AliAnalysisTaskSE(c),
 117   fBeamEnergy(c.fBeamEnergy),
 118   fInvMassFitLimitMin(c.fInvMassFitLimitMin),
 119   fInvMassFitLimitMax(c.fInvMassFitLimitMax),
 120   fPsiFitLimitMin(c.fPsiFitLimitMin),
 121   fPsiFitLimitMax(c.fPsiFitLimitMax),
 122   fPsiPFitLimitMin(c.fPsiPFitLimitMin),
 123   fPsiPFitLimitMax(c.fPsiPFitLimitMax),
 124   fBckFitLimitMin(c.fBckFitLimitMin),
 125   fBckFitLimitMax(c.fBckFitLimitMax),
 126   fInvariantMassFit(c.fInvariantMassFit),
 127   fkAnalysisType(c.fkAnalysisType),
 128   fOutput(c.fOutput)
 129  {
 130   //
 131   // Copy Constructor
 132   //
 133 }
 134
 135 //___________________________________________________________________________
 136 AliAnalysisTaskMuonDistributions::~AliAnalysisTaskMuonDistributions() {
 137   //
 138   //destructor
 139   //
 140   Info("~AliAnalysisTaskMuonDistributions","Calling Destructor");
 141
 142   if (fOutput){
 143     delete fOutput;
 144     fOutput = 0;
 145   }
 146 }
 147
 148 //___________________________________________________________________________
 149 void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
 150  //
 151  // output objects creation
 152  //
 153  fOutput = new TList();
 154  fOutput->SetOwner();
 155  //
 156  // various histos
 157  //
 158  TH1D *hNumberMuonTracks = new TH1D("hNumberMuonTracks","hNumberMuonTracks;N_{#mu tracks}",10,0.,10.);
 159  //
 160  // dimuon histos
 161  //
 162  TH1D *hMassDimu   = new TH1D("hMassDimu","hMassDimu;M_{#mu#mu} (GeV/c^{2})",180,0,9.);
 163  TH1D *hPtDimu  = new TH1D("hPtDimu","hPtDimu;p_{T} (GeV/c)",100,0,20);
 164  TH1D *hRapidityDimu  = new TH1D("hRapidityDimu","hRapidityDimu;y",100,-5,-2);
 165  TH1D *hCosThetaCSDimu  = new TH1D("hCosThetaCSDimu","hCosThetaCSDimu;cos#theta_{CS}",100,-1.,1.);
 166  TH1D *hCosThetaHEDimu  = new TH1D("hCosThetaHEDimu","hCosThetaHEDimu;cos#theta_{HE}",100,-1.,1.);
 167  //
 168  // muon histos
 169  //
 170  TH1D *hP  = new TH1D("hP","hP;p (GeV/c)",100,0,500);
 171  TH1D *hPt  = new TH1D("hPt","hPt;p_{T} (GeV/c)",100,0,20);
 172  TH1D *hRapidity  = new TH1D("hRapidity","hRapidity;y",100,-5,-2);
 173
 174  fOutput->Add(hNumberMuonTracks);
 175  fOutput->Add(hMassDimu);
 176  fOutput->Add(hPtDimu);
 177  fOutput->Add(hRapidityDimu);
 178  fOutput->Add(hCosThetaCSDimu);
 179  fOutput->Add(hCosThetaHEDimu);
 180  fOutput->Add(hP);
 181  fOutput->Add(hPt);
 182  fOutput->Add(hRapidity);
 183  fOutput->ls();
 184 }
 185
 186 //_________________________________________________
 187 void AliAnalysisTaskMuonDistributions::UserExec(Option_t *)
 188 {
 189 //
 190 // Execute analysis for current event
 191 //
 192   AliESDEvent *esd=0x0;
 193   AliAODEvent *aod=0x0;
 194
 195   if(strcmp(fkAnalysisType,"ESD")==0){
 196     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
 197         (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 198     if ( ! esdH ) {
 199       AliError("Cannot get input event handler");
 200       return;
 201     }
 202     esd = esdH->GetEvent();
 203   } else if(strcmp(fkAnalysisType,"AOD")==0){
 204     aod = dynamic_cast<AliAODEvent*> (InputEvent());
 205     if ( ! aod ) {
 206       AliError("Cannot get AOD event");
 207       return;
 208     }
 209   }
 210
 211   Int_t ntracks=-999;
 212   if(strcmp(fkAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks();
 213   else if(strcmp(fkAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks();
 214   Int_t nmuontracks=0;
 215
 216   for (Int_t j = 0; j<ntracks; j++) {
 217     Float_t pxmu1=-999; Float_t pymu1=-999; Float_t pzmu1=-999; Float_t ptmu1=-999; Float_t pmu1=-999;
 218     Float_t emu1=-999; Float_t rapiditymu1=-999;  Float_t chargemu1=-999;
 219     if(strcmp(fkAnalysisType,"ESD")==0){
 220       AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(esd->GetMuonTrack(j)));
 221       if (!mu1->ContainTrackerData()) continue;
 222       chargemu1 = mu1->Charge();
 223       pxmu1 = mu1->Px();
 224       pymu1 = mu1->Py();
 225       pzmu1 = mu1->Pz();
 226       emu1 = mu1->E();
 227       pmu1 = mu1->P();
 228       ptmu1 = mu1->Pt();
 229       rapiditymu1 = Rapidity(emu1,pzmu1);
 230     } else if(strcmp(fkAnalysisType,"AOD")==0){
 231       AliAODTrack *mu1 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
 232       if(!mu1) AliFatal("Not a standard AOD");
 233       if(!mu1->IsMuonTrack()) continue;
 234       chargemu1 = mu1->Charge();
 235       pxmu1 = mu1->Px();
 236       pymu1 = mu1->Py();
 237       pzmu1 = mu1->Pz();
 238       emu1 = mu1->E();
 239       pmu1 = mu1->P();
 240       ptmu1 = mu1->Pt();
 241       rapiditymu1 = Rapidity(emu1,pzmu1);
 242     }
 243     ((TH1D*)(fOutput->FindObject("hP")))->Fill(pmu1);
 244     ((TH1D*)(fOutput->FindObject("hPt")))->Fill(ptmu1);
 245     ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapiditymu1);
 246     nmuontracks++;
 247     if(chargemu1<0){
 248       for (Int_t jj = 0; jj<ntracks; jj++) {
 249         Float_t pxmu2=-999; Float_t pymu2=-999; Float_t pzmu2=-999;
 250         Float_t emu2=-999;Float_t chargemu2=-999;
 251         if(strcmp(fkAnalysisType,"ESD")==0){
 252           AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(esd->GetMuonTrack(jj)));
 253           if (!mu2->ContainTrackerData()) continue;
 254           chargemu2 = mu2->Charge();
 255           pxmu2 = mu2->Px();
 256           pymu2 = mu2->Py();
 257           pzmu2 = mu2->Pz();
 258           emu2 = mu2->E();
 259         } else if(strcmp(fkAnalysisType,"AOD")==0){
 260           AliAODTrack *mu2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(jj));
 261           if(!mu2) AliFatal("Not a standard AOD");
 262           if(!mu2->IsMuonTrack()) continue;
 263           chargemu2 = mu2->Charge();
 264           pxmu2 = mu2->Px();
 265           pymu2 = mu2->Py();
 266           pzmu2 = mu2->Pz();
 267           emu2 = mu2->E();
 268         }
 269         if(chargemu2>0){
 270           Float_t ptdimu = TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
 271           Float_t massdimu = InvMass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
 272           Float_t rapiditydimu = Rapidity((emu1+emu2),(pzmu1+pzmu2));
 273           Double_t costhetaCSdimu = CostCS((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
 274           Double_t costhetaHEdimu = CostHE((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
 275           ((TH1D*)(fOutput->FindObject("hMassDimu")))->Fill(massdimu);
 276           ((TH1D*)(fOutput->FindObject("hPtDimu")))->Fill(ptdimu);
 277           ((TH1D*)(fOutput->FindObject("hRapidityDimu")))->Fill(rapiditydimu);
 278           ((TH1D*)(fOutput->FindObject("hCosThetaCSDimu")))->Fill(costhetaCSdimu);
 279           ((TH1D*)(fOutput->FindObject("hCosThetaHEDimu")))->Fill(costhetaHEdimu);
 280         }
 281         //delete mu2;
 282       }      // second mu Loop
 283     }          // mu- Selection
 284     //delete mu1;
 285   }
 286   ((TH1D*)(fOutput->FindObject("hNumberMuonTracks")))->Fill(nmuontracks);
 287
 288   PostData(1,fOutput);
 289   }
 290
 291
 292 //________________________________________________________________________
 293 void AliAnalysisTaskMuonDistributions::Terminate(Option_t *)
 294 {
 295 //
 296 // Draw histos
 297 //
 298   gStyle->SetCanvasColor(10);
 299   gStyle->SetFrameFillColor(10);
 300   Int_t xmin=20;
 301   Int_t ymin=20;
 302
 303   printf("Using beam Energy=%f \n",fBeamEnergy);
 304
 305   fOutput = static_cast<TList*> (GetOutputData(1));
 306   TH1D *hNumberMuonTracks = static_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));
 307   TH1D *hMassDimu = static_cast<TH1D*> (fOutput->FindObject("hMassDimu"));
 308   TH1D *hPtDimu = static_cast<TH1D*> (fOutput->FindObject("hPtDimu"));
 309   TH1D *hRapidityDimu = static_cast<TH1D*> (fOutput->FindObject("hRapidityDimu"));
 310   TH1D *hCosThetaCSDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));
 311   TH1D *hCosThetaHEDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));
 312   TH1D *hP = static_cast<TH1D*> (fOutput->FindObject("hP"));
 313   TH1D *hPt = static_cast<TH1D*> (fOutput->FindObject("hPt"));
 314   TH1D *hRapidity = static_cast<TH1D*> (fOutput->FindObject("hRapidity"));
 315
 316
 317   TCanvas *c0_MuonDistributions = new TCanvas("c0_MuonDistributions","Plots",xmin,ymin,600,600);
 318   c0_MuonDistributions->Divide(2,2);
 319   c0_MuonDistributions->cd(1);
 320   hNumberMuonTracks->SetFillColor(2);
 321   hNumberMuonTracks->Draw();
 322
 323   xmin+=20; ymin+=20;
 324   TCanvas *c1_MuonDistributions = new TCanvas("c1_MuonDistributions","Muon kinematic distributions Plots",xmin,ymin,600,600);
 325   c1_MuonDistributions->Divide(2,2);
 326   c1_MuonDistributions->cd(1);
 327   gPad->SetLogy(1);
 328   hP->SetLineColor(4);
 329   hP->Draw();
 330   c1_MuonDistributions->cd(2);
 331   gPad->SetLogy(1);
 332   hPt->SetLineColor(4);
 333   hPt->Draw();
 334   c1_MuonDistributions->cd(3);
 335   hRapidity->SetLineColor(4);
 336   hRapidity->Draw();
 337
 338   xmin+=20; ymin+=20;
 339   TCanvas *c2_MuonDistributions = new TCanvas("c2_MuonDistributions","Dimuon kinematic distributions Plots",xmin,ymin,600,600);
 340   c2_MuonDistributions->Divide(2,2);
 341   c2_MuonDistributions->cd(1);
 342   gPad->SetLogy(1);
 343   hPtDimu->SetLineColor(2);
 344   hPtDimu->Draw();
 345   c2_MuonDistributions->cd(2);
 346   hRapidityDimu->SetLineColor(2);
 347   hRapidityDimu->Draw();
 348   c2_MuonDistributions->cd(3);
 349   hCosThetaCSDimu->SetLineColor(2);
 350   hCosThetaCSDimu->Draw();
 351   c2_MuonDistributions->cd(4);
 352   hCosThetaHEDimu->SetLineColor(2);
 353   hCosThetaHEDimu->Draw();
 354
 355   xmin+=20; ymin+=20;
 356   TCanvas *c3_MuonDistributions = new TCanvas("c3_MuonDistributions","Invariant Mass Plots",xmin,ymin,600,600);
 357   gPad->SetLogy(1);
 358   hMassDimu->Draw();
 359   if(fInvariantMassFit && hMassDimu->GetEntries()>100) FitInvMass(hMassDimu);
 360   c3_MuonDistributions->Update();
 361 }
 362
 363 //________________________________________________________________________
 364 Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 365                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 366 {
 367 //
 368 // invariant mass calculation
 369 //
 370     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 371                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 372     return imassrec;
 373 }
 374 //________________________________________________________________________
 375 Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const
 376 {
 377 //
 378 // calculate rapidity
 379 //
 380     Float_t rap;
 381     if(TMath::Abs(e-pz)>1e-7){
 382       rap = 0.5*TMath::Log((e+pz)/(e-pz));
 383       return rap;
 384     } else {
 385       rap = -200;
 386       return rap;
 387     }
 388 }
 389 //________________________________________________________________________
 390 Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 391 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 392 {
 393 //
 394 // costCS calculation
 395 //
 396   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 397   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 398   TVector3 beta,zaxisCS;
 399   Double_t mp=0.93827231;
 400   //
 401   // --- Fill the Lorentz vector for projectile and target in the CM frame
 402   //
 403   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 404   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 405   //
 406   // --- Get the muons parameters in the CM frame
 407   //
 408   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 409   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 410   //
 411   // --- Obtain the dimuon parameters in the CM frame
 412   //
 413   pDimuCM=pMu1CM+pMu2CM;
 414   //
 415   // --- Translate the dimuon parameters in the dimuon rest frame
 416   //
 417   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 418   pMu1Dimu=pMu1CM;
 419   pMu2Dimu=pMu2CM;
 420   pProjDimu=pProjCM;
 421   pTargDimu=pTargCM;
 422   pMu1Dimu.Boost(beta);
 423   pMu2Dimu.Boost(beta);
 424   pProjDimu.Boost(beta);
 425   pTargDimu.Boost(beta);
 426   //
 427   // --- Determine the z axis for the CS angle
 428   //
 429   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 430   //
 431   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 432   Double_t cost;
 433
 434   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 435   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 436   return cost;
 437 }
 438
 439 //________________________________________________________________________
 440 Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 441 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 442 {
 443 //
 444 // costHE calculation
 445 //
 446   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 447   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 448   TVector3 beta,zaxisCS;
 449   //
 450   // --- Get the muons parameters in the CM frame
 451   //
 452   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 453   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 454   //
 455   // --- Obtain the dimuon parameters in the CM frame
 456   //
 457   pDimuCM=pMu1CM+pMu2CM;
 458   //
 459   // --- Translate the muon parameters in the dimuon rest frame
 460   //
 461   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 462   pMu1Dimu=pMu1CM;
 463   pMu2Dimu=pMu2CM;
 464   pMu1Dimu.Boost(beta);
 465   pMu2Dimu.Boost(beta);
 466   //
 467   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 468   //
 469   TVector3 zaxis;
 470   zaxis=(pDimuCM.Vect()).Unit();
 471
 472   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 473   Double_t cost;
 474   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 475   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 476   return cost;
 477 }
 478
 479 //________________________________________________________________________
 480 void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
 481 {
 482 //
 483 // Fit to the Invariant Mass Spectrum
 484 //
 485   TF1 *gaupsi = new TF1("gaupsi","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
 486   TF1 *gaupsip = new TF1("gaupsip","gaus",fPsiPFitLimitMin,fPsiPFitLimitMax);
 487   TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);
 488   TF1 *tot = new TF1("mtot","gaus(0)+expo(3)+gaus(5)",fInvMassFitLimitMin,fInvMassFitLimitMax);
 489   Double_t par[8];
 490   Double_t binWidth= histo->GetBinWidth(1);
 491   gaupsi->SetLineColor(kGreen);
 492   gaupsi->SetLineWidth(2);
 493   histo->Fit(gaupsi,"Rl0");
 494   ex->SetLineColor(kBlue);
 495   ex->SetLineWidth(2);
 496   histo->Fit(ex,"Rl+");
 497   gaupsip->SetLineColor(kAzure-9);
 498   gaupsip->SetLineWidth(2);
 499   gaupsip->SetParLimits(1,3.6,3.75);
 500   gaupsip->SetParLimits(2,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2));
 501   histo->Fit(gaupsip,"Rl0+");
 502   gaupsi->GetParameters(&par[0]);
 503   ex->GetParameters(&par[3]);
 504   gaupsip->GetParameters(&par[5]);
 505   tot->SetParameters(par);
 506   tot->SetLineColor(2);
 507   tot->SetLineWidth(2);
 508   tot->SetParLimits(6,3.6,3.75);   //limit for the psi(2S) parameters
 509   tot->SetParLimits(7,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2)); //limit for the psi(2S) parameters
 510   histo->Fit(tot,"Rl+");
 511   histo->Draw("e");
 512   Double_t chi2 = tot->GetChisquare();
 513   Double_t ndf = tot->GetNDF();
 514
 515   Float_t meanPsi= tot->GetParameter(1);
 516   Float_t sigPsi= tot->GetParameter(2)*1000.;
 517   Double_t nPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/binWidth;
 518   TF1 *psifix = new TF1("psifix","gaus",2.,5.);
 519   psifix->SetParameter(0,tot->GetParameter(0));
 520   psifix->SetParameter(1,tot->GetParameter(1));
 521   psifix->SetParameter(2,tot->GetParameter(2));
 522   psifix->SetLineColor(kGreen);
 523   psifix->Draw("same");
 524   Double_t nPsi2933 = psifix->Integral(2.9,3.3)/binWidth;
 525
 526   TF1 *exfix = new TF1("exfix","expo",2.,5.);
 527   exfix->SetParameter(0,tot->GetParameter(3));
 528   exfix->SetParameter(1,tot->GetParameter(4));
 529   Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
 530
 531   Float_t meanPsiP= tot->GetParameter(6);
 532   Float_t sigPsiP= tot->GetParameter(7)*1000.;
 533   Double_t nPsiPFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(5)*tot->GetParameter(7)/binWidth;
 534   TF1 *psipfix = new TF1("psipfix","gaus",3.,5.);
 535   psipfix->SetParameter(0,tot->GetParameter(5));
 536   psipfix->SetParameter(1,tot->GetParameter(6));
 537   psipfix->SetParameter(2,tot->GetParameter(7));
 538   psipfix->SetLineColor(kAzure-9);
 539   psipfix->Draw("same");
 540
 541   printf("\n\n****************************************************************************\n");
 542   char psitext[100];
 543   if(nPsiFit<0) nPsiFit = 0;
 544   snprintf(psitext,100,"N. J/#psi = %10.0f",nPsiFit);
 545   printf("\nN. J/psi = %10.0f\n",nPsiFit);
 546   TLatex *psilatex = new TLatex(4.,0.85*histo->GetMaximum(),psitext);
 547   psilatex->SetTextColor(2);
 548   psilatex->SetTextSize(0.03);
 549   psilatex->SetTextAlign(2);
 550   psilatex->Draw();
 551
 552   char psi2text[100];
 553   snprintf(psi2text,100,"J/#psi m=%4.3f GeV   #sigma= %4.2f MeV",meanPsi,sigPsi);
 554   printf("J/psi m= %4.3f GeV sigma= %4.2f MeV\n",meanPsi,sigPsi);
 555   TLatex *psi2latex = new TLatex(4.,0.425*histo->GetMaximum(),psi2text);
 556   psi2latex->SetTextColor(2);
 557   psi2latex->SetTextSize(0.03);
 558   psi2latex->SetTextAlign(2);
 559   psi2latex->Draw();
 560
 561   char sbtext[100];
 562   snprintf(sbtext,100,"S/B (2.9-3.3)= %4.2f ",nPsi2933/nBck);
 563   printf("S/B (2.9-3.3)= %4.2f\n",nPsi2933/nBck);
 564   TLatex *sblatex = new TLatex(4.,0.212*histo->GetMaximum(),sbtext);
 565   sblatex->SetTextColor(2);
 566   sblatex->SetTextSize(0.03);
 567   sblatex->SetTextAlign(2);
 568   sblatex->Draw();
 569
 570   char psiptext[100];
 571   if(nPsiPFit<0) nPsiPFit = 0;
 572   snprintf(psiptext,100,"N. #psi(2S) = %10.0f",nPsiPFit);
 573   printf("\npsi(2S) = %10.0f\n",nPsiPFit);
 574   TLatex *psiplatex = new TLatex(4.,0.106*histo->GetMaximum(),psiptext);
 575   psiplatex->SetTextColor(2);
 576   psiplatex->SetTextSize(0.03);
 577   psiplatex->SetTextAlign(2);
 578   psiplatex->Draw();
 579
 580   char psip2text[100];
 581   snprintf(psip2text,100,"#psi(2S) m=%4.3f GeV   #sigma= %4.2f MeV",meanPsiP,sigPsiP);
 582   printf("psi(2S) m= %4.3f GeV sigma= %4.2f MeV\n",meanPsiP,sigPsiP);
 583   TLatex *psip2latex = new TLatex(4.,0.053*histo->GetMaximum(),psip2text);
 584   psip2latex->SetTextColor(2);
 585   psip2latex->SetTextSize(0.03);
 586   psip2latex->SetTextAlign(2);
 587   psip2latex->Draw();
 588
 589   char chi2text[100];
 590   snprintf(chi2text,100,"#chi^2/ndf = %4.2f ",chi2/ndf);
 591   printf("chi^2/ndf = %4.2f\n",chi2/ndf);
 592   TLatex *chi2latex = new TLatex(4.,0.026*histo->GetMaximum(),chi2text);
 593   chi2latex->SetTextColor(2);
 594   chi2latex->SetTextSize(0.03);
 595   chi2latex->SetTextAlign(2);
 596   chi2latex->Draw();
 597   printf("\n****************************************************************************\n");
 598
 599 }
 600
 601