]> 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
Minor fix for software triggers.
[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 <TCanvas.h>
96ee66b9 28#include <TF1.h>
e4dddba1 29#include <TH1.h>
e4dddba1 30#include <TLatex.h>
96ee66b9 31#include <TList.h>
32#include <TStyle.h>
e4dddba1 33
34#include "AliAnalysisTaskMuonDistributions.h"
35#include "AliAnalysisTaskSE.h"
36#include "AliAnalysisManager.h"
e4dddba1 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"
e4dddba1 51#include "AliESDInputHandler.h"
52#include "AliAODEvent.h"
53#include "AliAODHeader.h"
54#include "AliAODHandler.h"
55#include "AliAODInputHandler.h"
e4dddba1 56
57ClassImp(AliAnalysisTaskMuonDistributions)
58
59//__________________________________________________________________________
60AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions() :
61 fBeamEnergy(0.),
62 fInvMassFitLimitMin(2.),
63 fInvMassFitLimitMax(5.),
64 fPsiFitLimitMin(2.9),
65 fPsiFitLimitMax(3.3),
f032e1ac 66 fPsiPFitLimitMin(3.3),
67 fPsiPFitLimitMax(4.),
e4dddba1 68 fBckFitLimitMin(2.2),
69 fBckFitLimitMax(2.85),
70 fInvariantMassFit(kFALSE),
c5c5bf9c 71 fkAnalysisType(0x0),
e4dddba1 72 fOutput(0x0)
73{
74}
75//___________________________________________________________________________
76AliAnalysisTaskMuonDistributions::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),
f032e1ac 83 fPsiPFitLimitMin(3.3),
84 fPsiPFitLimitMax(4.),
e4dddba1 85 fBckFitLimitMin(2.2),
86 fBckFitLimitMax(2.85),
87 fInvariantMassFit(kFALSE),
c5c5bf9c 88 fkAnalysisType(0x0),
e4dddba1 89 fOutput(0x0)
90{
c5c5bf9c 91 //
e4dddba1 92 // Constructor. Initialization of Inputs and Outputs
93 //
94 Info("AliAnalysisTaskMuonDistributions","Calling Constructor");
95
96 DefineOutput(1,TList::Class());
97
98}
99
100//___________________________________________________________________________
101AliAnalysisTaskMuonDistributions& 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//___________________________________________________________________________
113AliAnalysisTaskMuonDistributions::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),
f032e1ac 120 fPsiPFitLimitMin(c.fPsiPFitLimitMin),
121 fPsiPFitLimitMax(c.fPsiPFitLimitMax),
e4dddba1 122 fBckFitLimitMin(c.fBckFitLimitMin),
123 fBckFitLimitMax(c.fBckFitLimitMax),
124 fInvariantMassFit(c.fInvariantMassFit),
c5c5bf9c 125 fkAnalysisType(c.fkAnalysisType),
e4dddba1 126 fOutput(c.fOutput)
127 {
128 //
129 // Copy Constructor
130 //
131}
132
133//___________________________________________________________________________
134AliAnalysisTaskMuonDistributions::~AliAnalysisTaskMuonDistributions() {
135 //
136 //destructor
137 //
138 Info("~AliAnalysisTaskMuonDistributions","Calling Destructor");
93fd3322 139
140 if (fOutput){
141 delete fOutput;
142 fOutput = 0;
143 }
e4dddba1 144}
145
146//___________________________________________________________________________
147void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
c5c5bf9c 148 //
149 // output objects creation
150 //
e4dddba1 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 //
c5c5bf9c 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.);
e4dddba1 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);
c5c5bf9c 173 fOutput->Add(hMassDimu);
174 fOutput->Add(hPtDimu);
175 fOutput->Add(hRapidityDimu);
176 fOutput->Add(hCosThetaCSDimu);
177 fOutput->Add(hCosThetaHEDimu);
e4dddba1 178 fOutput->Add(hP);
179 fOutput->Add(hPt);
180 fOutput->Add(hRapidity);
181 fOutput->ls();
182}
183
184//_________________________________________________
185void AliAnalysisTaskMuonDistributions::UserExec(Option_t *)
186{
c5c5bf9c 187//
188// Execute analysis for current event
189//
e4dddba1 190 AliESDEvent *esd=0x0;
191 AliAODEvent *aod=0x0;
192
c5c5bf9c 193 if(strcmp(fkAnalysisType,"ESD")==0){
e4dddba1 194 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
195 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
196 esd = esdH->GetEvent();
c5c5bf9c 197 } else if(strcmp(fkAnalysisType,"AOD")==0){
e4dddba1 198 aod = dynamic_cast<AliAODEvent*> (InputEvent());
199 }
200
201 Int_t ntracks=-999;
c5c5bf9c 202 if(strcmp(fkAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks();
203 else if(strcmp(fkAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks();
e4dddba1 204 Int_t nmuontracks=0;
205
206 for (Int_t j = 0; j<ntracks; j++) {
c5c5bf9c 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){
e4dddba1 210 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(esd->GetMuonTrack(j)));
211 if (!mu1->ContainTrackerData()) continue;
c5c5bf9c 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){
e4dddba1 221 AliAODTrack *mu1 = aod->GetTrack(j);
222 if(!mu1->IsMuonTrack()) continue;
c5c5bf9c 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);
e4dddba1 231 }
c5c5bf9c 232 ((TH1D*)(fOutput->FindObject("hP")))->Fill(pmu1);
233 ((TH1D*)(fOutput->FindObject("hPt")))->Fill(ptmu1);
234 ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapiditymu1);
e4dddba1 235 nmuontracks++;
c5c5bf9c 236 if(chargemu1<0){
e4dddba1 237 for (Int_t jj = 0; jj<ntracks; jj++) {
c5c5bf9c 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){
e4dddba1 241 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(esd->GetMuonTrack(jj)));
242 if (!mu2->ContainTrackerData()) continue;
c5c5bf9c 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){
e4dddba1 249 AliAODTrack *mu2 = aod->GetTrack(jj);
250 if(!mu2->IsMuonTrack()) continue;
c5c5bf9c 251 chargemu2 = mu2->Charge();
252 pxmu2 = mu2->Px();
253 pymu2 = mu2->Py();
254 pzmu2 = mu2->Pz();
255 emu2 = mu2->E();
e4dddba1 256 }
c5c5bf9c 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);
e4dddba1 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//________________________________________________________________________
281void AliAnalysisTaskMuonDistributions::Terminate(Option_t *)
282{
c5c5bf9c 283//
284// Draw histos
285//
e4dddba1 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
f032e1ac 293 fOutput = dynamic_cast<TList*> (GetOutputData(1));
c5c5bf9c 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"));
93fd3322 298 TH1D *hCosThetaCSDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));
299 TH1D *hCosThetaHEDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));
c5c5bf9c 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"));
f032e1ac 303
93fd3322 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();
524828c1 346 if(fInvariantMassFit && hMassDimu->GetEntries()>100) FitInvMass(hMassDimu);
93fd3322 347 c3_MuonDistributions->Update();
e4dddba1 348}
349
350//________________________________________________________________________
351Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
c5c5bf9c 352 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
e4dddba1 353{
c5c5bf9c 354//
e4dddba1 355// invariant mass calculation
c5c5bf9c 356//
e4dddba1 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//________________________________________________________________________
c5c5bf9c 362Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const
e4dddba1 363{
c5c5bf9c 364//
e4dddba1 365// calculate rapidity
c5c5bf9c 366//
e4dddba1 367 Float_t rap;
c5c5bf9c 368 if(TMath::Abs(e-pz)>1e-7){
e4dddba1 369 rap = 0.5*TMath::Log((e+pz)/(e-pz));
370 return rap;
371 } else {
372 rap = -200;
373 return rap;
374 }
375}
376//________________________________________________________________________
377Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
378Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
379{
f032e1ac 380//
c5c5bf9c 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
e4dddba1 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 //
c5c5bf9c 390 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
391 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
e4dddba1 392 //
393 // --- Get the muons parameters in the CM frame
394 //
c5c5bf9c 395 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
396 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
e4dddba1 397 //
398 // --- Obtain the dimuon parameters in the CM frame
399 //
c5c5bf9c 400 pDimuCM=pMu1CM+pMu2CM;
e4dddba1 401 //
402 // --- Translate the dimuon parameters in the dimuon rest frame
403 //
c5c5bf9c 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);
e4dddba1 413 //
414 // --- Determine the z axis for the CS angle
415 //
c5c5bf9c 416 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
e4dddba1 417 //
418 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
419 Double_t cost;
420
c5c5bf9c 421 if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
422 else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
e4dddba1 423 return cost;
424}
425
426//________________________________________________________________________
427Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
428Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
429{
c5c5bf9c 430//
431// costHE calculation
432//
433 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
434 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
e4dddba1 435 TVector3 beta,zaxisCS;
436 //
437 // --- Get the muons parameters in the CM frame
438 //
c5c5bf9c 439 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
440 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
e4dddba1 441 //
442 // --- Obtain the dimuon parameters in the CM frame
443 //
c5c5bf9c 444 pDimuCM=pMu1CM+pMu2CM;
e4dddba1 445 //
446 // --- Translate the muon parameters in the dimuon rest frame
447 //
c5c5bf9c 448 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
449 pMu1Dimu=pMu1CM;
450 pMu2Dimu=pMu2CM;
451 pMu1Dimu.Boost(beta);
452 pMu2Dimu.Boost(beta);
e4dddba1 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;
c5c5bf9c 457 zaxis=(pDimuCM.Vect()).Unit();
e4dddba1 458
459 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
460 Double_t cost;
c5c5bf9c 461 if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
462 else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
e4dddba1 463 return cost;
464}
465
466//________________________________________________________________________
467void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
468{
c5c5bf9c 469//
470// Fit to the Invariant Mass Spectrum
471//
f032e1ac 472 TF1 *gaupsi = new TF1("gaupsi","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
473 TF1 *gaupsip = new TF1("gaupsip","gaus",fPsiPFitLimitMin,fPsiPFitLimitMax);
e4dddba1 474 TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);
f032e1ac 475 TF1 *tot = new TF1("mtot","gaus(0)+expo(3)+gaus(5)",fInvMassFitLimitMin,fInvMassFitLimitMax);
476 Double_t par[8];
c5c5bf9c 477 Double_t binWidth= histo->GetBinWidth(1);
f032e1ac 478 gaupsi->SetLineColor(kGreen);
479 gaupsi->SetLineWidth(2);
480 histo->Fit(gaupsi,"Rl0");
481 ex->SetLineColor(kBlue);
e4dddba1 482 ex->SetLineWidth(2);
f032e1ac 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]);
e4dddba1 490 ex->GetParameters(&par[3]);
f032e1ac 491 gaupsip->GetParameters(&par[5]);
e4dddba1 492 tot->SetParameters(par);
493 tot->SetLineColor(2);
494 tot->SetLineWidth(2);
f032e1ac 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
e4dddba1 497 histo->Fit(tot,"Rl+");
498 histo->Draw("e");
c5c5bf9c 499 Double_t chi2 = tot->GetChisquare();
500 Double_t ndf = tot->GetNDF();
f032e1ac 501
c5c5bf9c 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;
f032e1ac 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
e4dddba1 513 TF1 *exfix = new TF1("exfix","expo",2.,5.);
514 exfix->SetParameter(0,tot->GetParameter(3));
515 exfix->SetParameter(1,tot->GetParameter(4));
c5c5bf9c 516 Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
e4dddba1 517
f032e1ac 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
e4dddba1 528 printf("\n\n****************************************************************************\n");
c5c5bf9c 529 char psitext[100];
524828c1 530 if(nPsiFit<0) nPsiFit = 0;
f032e1ac 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);
c5c5bf9c 534 psilatex->SetTextColor(2);
535 psilatex->SetTextSize(0.03);
536 psilatex->SetTextAlign(2);
537 psilatex->Draw();
e4dddba1 538
c5c5bf9c 539 char psi2text[100];
f032e1ac 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);
c5c5bf9c 543 psi2latex->SetTextColor(2);
544 psi2latex->SetTextSize(0.03);
545 psi2latex->SetTextAlign(2);
546 psi2latex->Draw();
e4dddba1 547
c5c5bf9c 548 char sbtext[100];
f032e1ac 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);
c5c5bf9c 552 sblatex->SetTextColor(2);
553 sblatex->SetTextSize(0.03);
554 sblatex->SetTextAlign(2);
555 sblatex->Draw();
e4dddba1 556
f032e1ac 557 char psiptext[100];
524828c1 558 if(nPsiPFit<0) nPsiPFit = 0;
f032e1ac 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
c5c5bf9c 576 char chi2text[100];
f032e1ac 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);
c5c5bf9c 580 chi2latex->SetTextColor(2);
581 chi2latex->SetTextSize(0.03);
582 chi2latex->SetTextAlign(2);
583 chi2latex->Draw();
e4dddba1 584 printf("\n****************************************************************************\n");
585
586}
587
588