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