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