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