]> 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
Changed AliAODEvent::GetTrack to return AliVTrack* (and not AliAODTrack)
[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){
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
2e8391b9 303 fOutput = static_cast<TList*> (GetOutputData(1));
304 TH1D *hNumberMuonTracks = static_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));
305 TH1D *hMassDimu = static_cast<TH1D*> (fOutput->FindObject("hMassDimu"));
306 TH1D *hPtDimu = static_cast<TH1D*> (fOutput->FindObject("hPtDimu"));
307 TH1D *hRapidityDimu = static_cast<TH1D*> (fOutput->FindObject("hRapidityDimu"));
308 TH1D *hCosThetaCSDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));
309 TH1D *hCosThetaHEDimu = static_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));
310 TH1D *hP = static_cast<TH1D*> (fOutput->FindObject("hP"));
311 TH1D *hPt = static_cast<TH1D*> (fOutput->FindObject("hPt"));
312 TH1D *hRapidity = static_cast<TH1D*> (fOutput->FindObject("hRapidity"));
313
f032e1ac 314
93fd3322 315 TCanvas *c0_MuonDistributions = new TCanvas("c0_MuonDistributions","Plots",xmin,ymin,600,600);
316 c0_MuonDistributions->Divide(2,2);
317 c0_MuonDistributions->cd(1);
318 hNumberMuonTracks->SetFillColor(2);
319 hNumberMuonTracks->Draw();
320
321 xmin+=20; ymin+=20;
322 TCanvas *c1_MuonDistributions = new TCanvas("c1_MuonDistributions","Muon kinematic distributions Plots",xmin,ymin,600,600);
323 c1_MuonDistributions->Divide(2,2);
324 c1_MuonDistributions->cd(1);
325 gPad->SetLogy(1);
326 hP->SetLineColor(4);
327 hP->Draw();
328 c1_MuonDistributions->cd(2);
329 gPad->SetLogy(1);
330 hPt->SetLineColor(4);
331 hPt->Draw();
332 c1_MuonDistributions->cd(3);
333 hRapidity->SetLineColor(4);
334 hRapidity->Draw();
335
336 xmin+=20; ymin+=20;
337 TCanvas *c2_MuonDistributions = new TCanvas("c2_MuonDistributions","Dimuon kinematic distributions Plots",xmin,ymin,600,600);
338 c2_MuonDistributions->Divide(2,2);
339 c2_MuonDistributions->cd(1);
340 gPad->SetLogy(1);
341 hPtDimu->SetLineColor(2);
342 hPtDimu->Draw();
343 c2_MuonDistributions->cd(2);
344 hRapidityDimu->SetLineColor(2);
345 hRapidityDimu->Draw();
346 c2_MuonDistributions->cd(3);
347 hCosThetaCSDimu->SetLineColor(2);
348 hCosThetaCSDimu->Draw();
349 c2_MuonDistributions->cd(4);
350 hCosThetaHEDimu->SetLineColor(2);
351 hCosThetaHEDimu->Draw();
352
353 xmin+=20; ymin+=20;
354 TCanvas *c3_MuonDistributions = new TCanvas("c3_MuonDistributions","Invariant Mass Plots",xmin,ymin,600,600);
355 gPad->SetLogy(1);
356 hMassDimu->Draw();
524828c1 357 if(fInvariantMassFit && hMassDimu->GetEntries()>100) FitInvMass(hMassDimu);
93fd3322 358 c3_MuonDistributions->Update();
e4dddba1 359}
360
361//________________________________________________________________________
362Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
c5c5bf9c 363 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
e4dddba1 364{
c5c5bf9c 365//
e4dddba1 366// invariant mass calculation
c5c5bf9c 367//
e4dddba1 368 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
369 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
370 return imassrec;
371}
372//________________________________________________________________________
c5c5bf9c 373Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const
e4dddba1 374{
c5c5bf9c 375//
e4dddba1 376// calculate rapidity
c5c5bf9c 377//
e4dddba1 378 Float_t rap;
c5c5bf9c 379 if(TMath::Abs(e-pz)>1e-7){
e4dddba1 380 rap = 0.5*TMath::Log((e+pz)/(e-pz));
381 return rap;
382 } else {
383 rap = -200;
384 return rap;
385 }
386}
387//________________________________________________________________________
388Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
389Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
390{
f032e1ac 391//
c5c5bf9c 392// costCS calculation
393//
394 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
395 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
e4dddba1 396 TVector3 beta,zaxisCS;
397 Double_t mp=0.93827231;
398 //
399 // --- Fill the Lorentz vector for projectile and target in the CM frame
400 //
c5c5bf9c 401 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
402 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
e4dddba1 403 //
404 // --- Get the muons parameters in the CM frame
405 //
c5c5bf9c 406 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
407 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
e4dddba1 408 //
409 // --- Obtain the dimuon parameters in the CM frame
410 //
c5c5bf9c 411 pDimuCM=pMu1CM+pMu2CM;
e4dddba1 412 //
413 // --- Translate the dimuon parameters in the dimuon rest frame
414 //
c5c5bf9c 415 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
416 pMu1Dimu=pMu1CM;
417 pMu2Dimu=pMu2CM;
418 pProjDimu=pProjCM;
419 pTargDimu=pTargCM;
420 pMu1Dimu.Boost(beta);
421 pMu2Dimu.Boost(beta);
422 pProjDimu.Boost(beta);
423 pTargDimu.Boost(beta);
e4dddba1 424 //
425 // --- Determine the z axis for the CS angle
426 //
c5c5bf9c 427 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
e4dddba1 428 //
429 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
430 Double_t cost;
431
c5c5bf9c 432 if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
433 else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
e4dddba1 434 return cost;
435}
436
437//________________________________________________________________________
438Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
439Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
440{
c5c5bf9c 441//
442// costHE calculation
443//
444 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
445 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
e4dddba1 446 TVector3 beta,zaxisCS;
447 //
448 // --- Get the muons parameters in the CM frame
449 //
c5c5bf9c 450 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
451 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
e4dddba1 452 //
453 // --- Obtain the dimuon parameters in the CM frame
454 //
c5c5bf9c 455 pDimuCM=pMu1CM+pMu2CM;
e4dddba1 456 //
457 // --- Translate the muon parameters in the dimuon rest frame
458 //
c5c5bf9c 459 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
460 pMu1Dimu=pMu1CM;
461 pMu2Dimu=pMu2CM;
462 pMu1Dimu.Boost(beta);
463 pMu2Dimu.Boost(beta);
e4dddba1 464 //
465 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
466 //
467 TVector3 zaxis;
c5c5bf9c 468 zaxis=(pDimuCM.Vect()).Unit();
e4dddba1 469
470 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
471 Double_t cost;
c5c5bf9c 472 if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
473 else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
e4dddba1 474 return cost;
475}
476
477//________________________________________________________________________
478void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
479{
c5c5bf9c 480//
481// Fit to the Invariant Mass Spectrum
482//
f032e1ac 483 TF1 *gaupsi = new TF1("gaupsi","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
484 TF1 *gaupsip = new TF1("gaupsip","gaus",fPsiPFitLimitMin,fPsiPFitLimitMax);
e4dddba1 485 TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);
f032e1ac 486 TF1 *tot = new TF1("mtot","gaus(0)+expo(3)+gaus(5)",fInvMassFitLimitMin,fInvMassFitLimitMax);
487 Double_t par[8];
c5c5bf9c 488 Double_t binWidth= histo->GetBinWidth(1);
f032e1ac 489 gaupsi->SetLineColor(kGreen);
490 gaupsi->SetLineWidth(2);
491 histo->Fit(gaupsi,"Rl0");
492 ex->SetLineColor(kBlue);
e4dddba1 493 ex->SetLineWidth(2);
f032e1ac 494 histo->Fit(ex,"Rl+");
495 gaupsip->SetLineColor(kAzure-9);
496 gaupsip->SetLineWidth(2);
497 gaupsip->SetParLimits(1,3.6,3.75);
498 gaupsip->SetParLimits(2,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2));
499 histo->Fit(gaupsip,"Rl0+");
500 gaupsi->GetParameters(&par[0]);
e4dddba1 501 ex->GetParameters(&par[3]);
f032e1ac 502 gaupsip->GetParameters(&par[5]);
e4dddba1 503 tot->SetParameters(par);
504 tot->SetLineColor(2);
505 tot->SetLineWidth(2);
f032e1ac 506 tot->SetParLimits(6,3.6,3.75); //limit for the psi(2S) parameters
507 tot->SetParLimits(7,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2)); //limit for the psi(2S) parameters
e4dddba1 508 histo->Fit(tot,"Rl+");
509 histo->Draw("e");
c5c5bf9c 510 Double_t chi2 = tot->GetChisquare();
511 Double_t ndf = tot->GetNDF();
f032e1ac 512
c5c5bf9c 513 Float_t meanPsi= tot->GetParameter(1);
514 Float_t sigPsi= tot->GetParameter(2)*1000.;
515 Double_t nPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/binWidth;
f032e1ac 516 TF1 *psifix = new TF1("psifix","gaus",2.,5.);
517 psifix->SetParameter(0,tot->GetParameter(0));
518 psifix->SetParameter(1,tot->GetParameter(1));
519 psifix->SetParameter(2,tot->GetParameter(2));
520 psifix->SetLineColor(kGreen);
521 psifix->Draw("same");
522 Double_t nPsi2933 = psifix->Integral(2.9,3.3)/binWidth;
523
e4dddba1 524 TF1 *exfix = new TF1("exfix","expo",2.,5.);
525 exfix->SetParameter(0,tot->GetParameter(3));
526 exfix->SetParameter(1,tot->GetParameter(4));
c5c5bf9c 527 Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
e4dddba1 528
f032e1ac 529 Float_t meanPsiP= tot->GetParameter(6);
530 Float_t sigPsiP= tot->GetParameter(7)*1000.;
531 Double_t nPsiPFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(5)*tot->GetParameter(7)/binWidth;
532 TF1 *psipfix = new TF1("psipfix","gaus",3.,5.);
533 psipfix->SetParameter(0,tot->GetParameter(5));
534 psipfix->SetParameter(1,tot->GetParameter(6));
535 psipfix->SetParameter(2,tot->GetParameter(7));
536 psipfix->SetLineColor(kAzure-9);
537 psipfix->Draw("same");
538
e4dddba1 539 printf("\n\n****************************************************************************\n");
c5c5bf9c 540 char psitext[100];
524828c1 541 if(nPsiFit<0) nPsiFit = 0;
7bd97ad1 542 snprintf(psitext,100,"N. J/#psi = %10.0f",nPsiFit);
f032e1ac 543 printf("\nN. J/psi = %10.0f\n",nPsiFit);
544 TLatex *psilatex = new TLatex(4.,0.85*histo->GetMaximum(),psitext);
c5c5bf9c 545 psilatex->SetTextColor(2);
546 psilatex->SetTextSize(0.03);
547 psilatex->SetTextAlign(2);
548 psilatex->Draw();
e4dddba1 549
c5c5bf9c 550 char psi2text[100];
7bd97ad1 551 snprintf(psi2text,100,"J/#psi m=%4.3f GeV #sigma= %4.2f MeV",meanPsi,sigPsi);
f032e1ac 552 printf("J/psi m= %4.3f GeV sigma= %4.2f MeV\n",meanPsi,sigPsi);
553 TLatex *psi2latex = new TLatex(4.,0.425*histo->GetMaximum(),psi2text);
c5c5bf9c 554 psi2latex->SetTextColor(2);
555 psi2latex->SetTextSize(0.03);
556 psi2latex->SetTextAlign(2);
557 psi2latex->Draw();
e4dddba1 558
c5c5bf9c 559 char sbtext[100];
7bd97ad1 560 snprintf(sbtext,100,"S/B (2.9-3.3)= %4.2f ",nPsi2933/nBck);
f032e1ac 561 printf("S/B (2.9-3.3)= %4.2f\n",nPsi2933/nBck);
562 TLatex *sblatex = new TLatex(4.,0.212*histo->GetMaximum(),sbtext);
c5c5bf9c 563 sblatex->SetTextColor(2);
564 sblatex->SetTextSize(0.03);
565 sblatex->SetTextAlign(2);
566 sblatex->Draw();
e4dddba1 567
f032e1ac 568 char psiptext[100];
524828c1 569 if(nPsiPFit<0) nPsiPFit = 0;
7bd97ad1 570 snprintf(psiptext,100,"N. #psi(2S) = %10.0f",nPsiPFit);
f032e1ac 571 printf("\npsi(2S) = %10.0f\n",nPsiPFit);
572 TLatex *psiplatex = new TLatex(4.,0.106*histo->GetMaximum(),psiptext);
573 psiplatex->SetTextColor(2);
574 psiplatex->SetTextSize(0.03);
575 psiplatex->SetTextAlign(2);
576 psiplatex->Draw();
577
578 char psip2text[100];
7bd97ad1 579 snprintf(psip2text,100,"#psi(2S) m=%4.3f GeV #sigma= %4.2f MeV",meanPsiP,sigPsiP);
f032e1ac 580 printf("psi(2S) m= %4.3f GeV sigma= %4.2f MeV\n",meanPsiP,sigPsiP);
581 TLatex *psip2latex = new TLatex(4.,0.053*histo->GetMaximum(),psip2text);
582 psip2latex->SetTextColor(2);
583 psip2latex->SetTextSize(0.03);
584 psip2latex->SetTextAlign(2);
585 psip2latex->Draw();
586
c5c5bf9c 587 char chi2text[100];
7bd97ad1 588 snprintf(chi2text,100,"#chi^2/ndf = %4.2f ",chi2/ndf);
f032e1ac 589 printf("chi^2/ndf = %4.2f\n",chi2/ndf);
590 TLatex *chi2latex = new TLatex(4.,0.026*histo->GetMaximum(),chi2text);
c5c5bf9c 591 chi2latex->SetTextColor(2);
592 chi2latex->SetTextSize(0.03);
593 chi2latex->SetTextAlign(2);
594 chi2latex->Draw();
e4dddba1 595 printf("\n****************************************************************************\n");
596
597}
598
599