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 = aod->GetTrack(j);
 232       if(!mu1->IsMuonTrack()) continue;
 233       chargemu1 = mu1->Charge();
 234       pxmu1 = mu1->Px();
 235       pymu1 = mu1->Py();
 236       pzmu1 = mu1->Pz();
 237       emu1 = mu1->E();
 238       pmu1 = mu1->P();
 239       ptmu1 = mu1->Pt();
 240       rapiditymu1 = Rapidity(emu1,pzmu1);
 241     }
 242     ((TH1D*)(fOutput->FindObject("hP")))->Fill(pmu1);
 243     ((TH1D*)(fOutput->FindObject("hPt")))->Fill(ptmu1);
 244     ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapiditymu1);
 245     nmuontracks++;
 246     if(chargemu1<0){
 247       for (Int_t jj = 0; jj<ntracks; jj++) {
 248         Float_t pxmu2=-999; Float_t pymu2=-999; Float_t pzmu2=-999;
 249         Float_t emu2=-999;Float_t chargemu2=-999;
 250         if(strcmp(fkAnalysisType,"ESD")==0){
 251           AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(esd->GetMuonTrack(jj)));
 252           if (!mu2->ContainTrackerData()) continue;
 253           chargemu2 = mu2->Charge();
 254           pxmu2 = mu2->Px();
 255           pymu2 = mu2->Py();
 256           pzmu2 = mu2->Pz();
 257           emu2 = mu2->E();
 258         } else if(strcmp(fkAnalysisType,"AOD")==0){
 259           AliAODTrack *mu2 = aod->GetTrack(jj);
 260           if(!mu2->IsMuonTrack()) continue;
 261           chargemu2 = mu2->Charge();
 262           pxmu2 = mu2->Px();
 263           pymu2 = mu2->Py();
 264           pzmu2 = mu2->Pz();
 265           emu2 = mu2->E();
 266         }
 267         if(chargemu2>0){
 268           Float_t ptdimu = TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
 269           Float_t massdimu = InvMass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
 270           Float_t rapiditydimu = Rapidity((emu1+emu2),(pzmu1+pzmu2));
 271           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);
 272           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);
 273           ((TH1D*)(fOutput->FindObject("hMassDimu")))->Fill(massdimu);
 274           ((TH1D*)(fOutput->FindObject("hPtDimu")))->Fill(ptdimu);
 275           ((TH1D*)(fOutput->FindObject("hRapidityDimu")))->Fill(rapiditydimu);
 276           ((TH1D*)(fOutput->FindObject("hCosThetaCSDimu")))->Fill(costhetaCSdimu);
 277           ((TH1D*)(fOutput->FindObject("hCosThetaHEDimu")))->Fill(costhetaHEdimu);
 278         }
 279         //delete mu2;
 280       }      // second mu Loop
 281     }          // mu- Selection
 282     //delete mu1;
 283   }
 284   ((TH1D*)(fOutput->FindObject("hNumberMuonTracks")))->Fill(nmuontracks);
 285
 286   PostData(1,fOutput);
 287   }
 288
 289
 290 //________________________________________________________________________
 291 void AliAnalysisTaskMuonDistributions::Terminate(Option_t *)
 292 {
 293 //
 294 // Draw histos
 295 //
 296   gStyle->SetCanvasColor(10);
 297   gStyle->SetFrameFillColor(10);
 298   Int_t xmin=20;
 299   Int_t ymin=20;
 300
 301   printf("Using beam Energy=%f \n",fBeamEnergy);
 302
 303   fOutput = static_cast<TList*> (GetOutputData(1));
 304   TH1D *hNumberMuonTracks = static_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));
 305   TH1D *hMassDimu = static_cast<TH1D*> (fOutput->FindObject("hMassDimu"));
 306   TH1D *hPtDimu = static_cast<TH1D*> (fOutput->FindObject("hPtDimu"));
 307   TH1D *hRapidityDimu = static_cast<TH1D*> (fOutput->FindObject("hRapidityDimu"));
 308   TH1D *hCosThetaCSDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));
 309   TH1D *hCosThetaHEDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));
 310   TH1D *hP = static_cast<TH1D*> (fOutput->FindObject("hP"));
 311   TH1D *hPt = static_cast<TH1D*> (fOutput->FindObject("hPt"));
 312   TH1D *hRapidity = static_cast<TH1D*> (fOutput->FindObject("hRapidity"));
 313
 314
 315   TCanvas *c0_MuonDistributions = new TCanvas("c0_MuonDistributions","Plots",xmin,ymin,600,600);
 316   c0_MuonDistributions->Divide(2,2);
 317   c0_MuonDistributions->cd(1);
 318   hNumberMuonTracks->SetFillColor(2);
 319   hNumberMuonTracks->Draw();
 320
 321   xmin+=20; ymin+=20;
 322   TCanvas *c1_MuonDistributions = new TCanvas("c1_MuonDistributions","Muon kinematic distributions Plots",xmin,ymin,600,600);
 323   c1_MuonDistributions->Divide(2,2);
 324   c1_MuonDistributions->cd(1);
 325   gPad->SetLogy(1);
 326   hP->SetLineColor(4);
 327   hP->Draw();
 328   c1_MuonDistributions->cd(2);
 329   gPad->SetLogy(1);
 330   hPt->SetLineColor(4);
 331   hPt->Draw();
 332   c1_MuonDistributions->cd(3);
 333   hRapidity->SetLineColor(4);
 334   hRapidity->Draw();
 335
 336   xmin+=20; ymin+=20;
 337   TCanvas *c2_MuonDistributions = new TCanvas("c2_MuonDistributions","Dimuon kinematic distributions Plots",xmin,ymin,600,600);
 338   c2_MuonDistributions->Divide(2,2);
 339   c2_MuonDistributions->cd(1);
 340   gPad->SetLogy(1);
 341   hPtDimu->SetLineColor(2);
 342   hPtDimu->Draw();
 343   c2_MuonDistributions->cd(2);
 344   hRapidityDimu->SetLineColor(2);
 345   hRapidityDimu->Draw();
 346   c2_MuonDistributions->cd(3);
 347   hCosThetaCSDimu->SetLineColor(2);
 348   hCosThetaCSDimu->Draw();
 349   c2_MuonDistributions->cd(4);
 350   hCosThetaHEDimu->SetLineColor(2);
 351   hCosThetaHEDimu->Draw();
 352
 353   xmin+=20; ymin+=20;
 354   TCanvas *c3_MuonDistributions = new TCanvas("c3_MuonDistributions","Invariant Mass Plots",xmin,ymin,600,600);
 355   gPad->SetLogy(1);
 356   hMassDimu->Draw();
 357   if(fInvariantMassFit && hMassDimu->GetEntries()>100) FitInvMass(hMassDimu);
 358   c3_MuonDistributions->Update();
 359 }
 360
 361 //________________________________________________________________________
 362 Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 363                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 364 {
 365 //
 366 // invariant mass calculation
 367 //
 368     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 369                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 370     return imassrec;
 371 }
 372 //________________________________________________________________________
 373 Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const
 374 {
 375 //
 376 // calculate rapidity
 377 //
 378     Float_t rap;
 379     if(TMath::Abs(e-pz)>1e-7){
 380       rap = 0.5*TMath::Log((e+pz)/(e-pz));
 381       return rap;
 382     } else {
 383       rap = -200;
 384       return rap;
 385     }
 386 }
 387 //________________________________________________________________________
 388 Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 389 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 390 {
 391 //
 392 // costCS calculation
 393 //
 394   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 395   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 396   TVector3 beta,zaxisCS;
 397   Double_t mp=0.93827231;
 398   //
 399   // --- Fill the Lorentz vector for projectile and target in the CM frame
 400   //
 401   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 402   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 403   //
 404   // --- Get the muons parameters in the CM frame
 405   //
 406   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 407   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 408   //
 409   // --- Obtain the dimuon parameters in the CM frame
 410   //
 411   pDimuCM=pMu1CM+pMu2CM;
 412   //
 413   // --- Translate the dimuon parameters in the dimuon rest frame
 414   //
 415   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 416   pMu1Dimu=pMu1CM;
 417   pMu2Dimu=pMu2CM;
 418   pProjDimu=pProjCM;
 419   pTargDimu=pTargCM;
 420   pMu1Dimu.Boost(beta);
 421   pMu2Dimu.Boost(beta);
 422   pProjDimu.Boost(beta);
 423   pTargDimu.Boost(beta);
 424   //
 425   // --- Determine the z axis for the CS angle
 426   //
 427   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 428   //
 429   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 430   Double_t cost;
 431
 432   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 433   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 434   return cost;
 435 }
 436
 437 //________________________________________________________________________
 438 Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 439 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 440 {
 441 //
 442 // costHE calculation
 443 //
 444   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 445   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 446   TVector3 beta,zaxisCS;
 447   //
 448   // --- Get the muons parameters in the CM frame
 449   //
 450   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 451   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 452   //
 453   // --- Obtain the dimuon parameters in the CM frame
 454   //
 455   pDimuCM=pMu1CM+pMu2CM;
 456   //
 457   // --- Translate the muon parameters in the dimuon rest frame
 458   //
 459   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 460   pMu1Dimu=pMu1CM;
 461   pMu2Dimu=pMu2CM;
 462   pMu1Dimu.Boost(beta);
 463   pMu2Dimu.Boost(beta);
 464   //
 465   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 466   //
 467   TVector3 zaxis;
 468   zaxis=(pDimuCM.Vect()).Unit();
 469
 470   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 471   Double_t cost;
 472   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 473   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 474   return cost;
 475 }
 476
 477 //________________________________________________________________________
 478 void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
 479 {
 480 //
 481 // Fit to the Invariant Mass Spectrum
 482 //
 483   TF1 *gaupsi = new TF1("gaupsi","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
 484   TF1 *gaupsip = new TF1("gaupsip","gaus",fPsiPFitLimitMin,fPsiPFitLimitMax);
 485   TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);
 486   TF1 *tot = new TF1("mtot","gaus(0)+expo(3)+gaus(5)",fInvMassFitLimitMin,fInvMassFitLimitMax);
 487   Double_t par[8];
 488   Double_t binWidth= histo->GetBinWidth(1);
 489   gaupsi->SetLineColor(kGreen);
 490   gaupsi->SetLineWidth(2);
 491   histo->Fit(gaupsi,"Rl0");
 492   ex->SetLineColor(kBlue);
 493   ex->SetLineWidth(2);
 494   histo->Fit(ex,"Rl+");
 495   gaupsip->SetLineColor(kAzure-9);
 496   gaupsip->SetLineWidth(2);
 497   gaupsip->SetParLimits(1,3.6,3.75);
 498   gaupsip->SetParLimits(2,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2));
 499   histo->Fit(gaupsip,"Rl0+");
 500   gaupsi->GetParameters(&par[0]);
 501   ex->GetParameters(&par[3]);
 502   gaupsip->GetParameters(&par[5]);
 503   tot->SetParameters(par);
 504   tot->SetLineColor(2);
 505   tot->SetLineWidth(2);
 506   tot->SetParLimits(6,3.6,3.75);   //limit for the psi(2S) parameters
 507   tot->SetParLimits(7,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2)); //limit for the psi(2S) parameters
 508   histo->Fit(tot,"Rl+");
 509   histo->Draw("e");
 510   Double_t chi2 = tot->GetChisquare();
 511   Double_t ndf = tot->GetNDF();
 512
 513   Float_t meanPsi= tot->GetParameter(1);
 514   Float_t sigPsi= tot->GetParameter(2)*1000.;
 515   Double_t nPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/binWidth;
 516   TF1 *psifix = new TF1("psifix","gaus",2.,5.);
 517   psifix->SetParameter(0,tot->GetParameter(0));
 518   psifix->SetParameter(1,tot->GetParameter(1));
 519   psifix->SetParameter(2,tot->GetParameter(2));
 520   psifix->SetLineColor(kGreen);
 521   psifix->Draw("same");
 522   Double_t nPsi2933 = psifix->Integral(2.9,3.3)/binWidth;
 523
 524   TF1 *exfix = new TF1("exfix","expo",2.,5.);
 525   exfix->SetParameter(0,tot->GetParameter(3));
 526   exfix->SetParameter(1,tot->GetParameter(4));
 527   Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
 528
 529   Float_t meanPsiP= tot->GetParameter(6);
 530   Float_t sigPsiP= tot->GetParameter(7)*1000.;
 531   Double_t nPsiPFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(5)*tot->GetParameter(7)/binWidth;
 532   TF1 *psipfix = new TF1("psipfix","gaus",3.,5.);
 533   psipfix->SetParameter(0,tot->GetParameter(5));
 534   psipfix->SetParameter(1,tot->GetParameter(6));
 535   psipfix->SetParameter(2,tot->GetParameter(7));
 536   psipfix->SetLineColor(kAzure-9);
 537   psipfix->Draw("same");
 538
 539   printf("\n\n****************************************************************************\n");
 540   char psitext[100];
 541   if(nPsiFit<0) nPsiFit = 0;
 542   snprintf(psitext,100,"N. J/#psi = %10.0f",nPsiFit);
 543   printf("\nN. J/psi = %10.0f\n",nPsiFit);
 544   TLatex *psilatex = new TLatex(4.,0.85*histo->GetMaximum(),psitext);
 545   psilatex->SetTextColor(2);
 546   psilatex->SetTextSize(0.03);
 547   psilatex->SetTextAlign(2);
 548   psilatex->Draw();
 549
 550   char psi2text[100];
 551   snprintf(psi2text,100,"J/#psi m=%4.3f GeV   #sigma= %4.2f MeV",meanPsi,sigPsi);
 552   printf("J/psi m= %4.3f GeV sigma= %4.2f MeV\n",meanPsi,sigPsi);
 553   TLatex *psi2latex = new TLatex(4.,0.425*histo->GetMaximum(),psi2text);
 554   psi2latex->SetTextColor(2);
 555   psi2latex->SetTextSize(0.03);
 556   psi2latex->SetTextAlign(2);
 557   psi2latex->Draw();
 558
 559   char sbtext[100];
 560   snprintf(sbtext,100,"S/B (2.9-3.3)= %4.2f ",nPsi2933/nBck);
 561   printf("S/B (2.9-3.3)= %4.2f\n",nPsi2933/nBck);
 562   TLatex *sblatex = new TLatex(4.,0.212*histo->GetMaximum(),sbtext);
 563   sblatex->SetTextColor(2);
 564   sblatex->SetTextSize(0.03);
 565   sblatex->SetTextAlign(2);
 566   sblatex->Draw();
 567
 568   char psiptext[100];
 569   if(nPsiPFit<0) nPsiPFit = 0;
 570   snprintf(psiptext,100,"N. #psi(2S) = %10.0f",nPsiPFit);
 571   printf("\npsi(2S) = %10.0f\n",nPsiPFit);
 572   TLatex *psiplatex = new TLatex(4.,0.106*histo->GetMaximum(),psiptext);
 573   psiplatex->SetTextColor(2);
 574   psiplatex->SetTextSize(0.03);
 575   psiplatex->SetTextAlign(2);
 576   psiplatex->Draw();
 577
 578   char psip2text[100];
 579   snprintf(psip2text,100,"#psi(2S) m=%4.3f GeV   #sigma= %4.2f MeV",meanPsiP,sigPsiP);
 580   printf("psi(2S) m= %4.3f GeV sigma= %4.2f MeV\n",meanPsiP,sigPsiP);
 581   TLatex *psip2latex = new TLatex(4.,0.053*histo->GetMaximum(),psip2text);
 582   psip2latex->SetTextColor(2);
 583   psip2latex->SetTextSize(0.03);
 584   psip2latex->SetTextAlign(2);
 585   psip2latex->Draw();
 586
 587   char chi2text[100];
 588   snprintf(chi2text,100,"#chi^2/ndf = %4.2f ",chi2/ndf);
 589   printf("chi^2/ndf = %4.2f\n",chi2/ndf);
 590   TLatex *chi2latex = new TLatex(4.,0.026*histo->GetMaximum(),chi2text);
 591   chi2latex->SetTextColor(2);
 592   chi2latex->SetTextSize(0.03);
 593   chi2latex->SetTextAlign(2);
 594   chi2latex->Draw();
 595   printf("\n****************************************************************************\n");
 596
 597 }
 598
 599