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