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