1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------------------
17 // Analysis task to compute muon/dimuon kinematic distributions
18 // The output is a list of histograms.
19 // The macro class can run on AOD or in the train with the ESD filter.
22 //-----------------------------------------------------------------------------
24 //#ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
25 //#define ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
34 #include "AliAnalysisTaskMuonDistributions.h"
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisManager.h"
37 #include "AliESDEvent.h"
39 #include "AliVEvent.h"
40 #include "AliMCEventHandler.h"
41 #include "AliInputEventHandler.h"
42 #include "AliMCEvent.h"
45 #include "AliHeader.h"
46 #include "AliESDHeader.h"
48 #include "TParticle.h"
49 #include "TLorentzVector.h"
50 #include "AliESDMuonTrack.h"
51 #include "AliESDInputHandler.h"
52 #include "AliAODEvent.h"
53 #include "AliAODHeader.h"
54 #include "AliAODHandler.h"
55 #include "AliAODInputHandler.h"
57 ClassImp(AliAnalysisTaskMuonDistributions)
59 //__________________________________________________________________________
60 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions() :
62 fInvMassFitLimitMin(2.),
63 fInvMassFitLimitMax(5.),
66 fPsiPFitLimitMin(3.3),
69 fBckFitLimitMax(2.85),
70 fInvariantMassFit(kFALSE),
75 //___________________________________________________________________________
76 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const Char_t* name) :
77 AliAnalysisTaskSE(name),
79 fInvMassFitLimitMin(2.),
80 fInvMassFitLimitMax(5.),
83 fPsiPFitLimitMin(3.3),
86 fBckFitLimitMax(2.85),
87 fInvariantMassFit(kFALSE),
92 // Constructor. Initialization of Inputs and Outputs
94 Info("AliAnalysisTaskMuonDistributions","Calling Constructor");
96 DefineOutput(1,TList::Class());
100 //___________________________________________________________________________
101 AliAnalysisTaskMuonDistributions& AliAnalysisTaskMuonDistributions::operator=(const AliAnalysisTaskMuonDistributions& c)
104 // Assignment operator
107 AliAnalysisTaskSE::operator=(c) ;
112 //___________________________________________________________________________
113 AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const AliAnalysisTaskMuonDistributions& c) :
114 AliAnalysisTaskSE(c),
115 fBeamEnergy(c.fBeamEnergy),
116 fInvMassFitLimitMin(c.fInvMassFitLimitMin),
117 fInvMassFitLimitMax(c.fInvMassFitLimitMax),
118 fPsiFitLimitMin(c.fPsiFitLimitMin),
119 fPsiFitLimitMax(c.fPsiFitLimitMax),
120 fPsiPFitLimitMin(c.fPsiPFitLimitMin),
121 fPsiPFitLimitMax(c.fPsiPFitLimitMax),
122 fBckFitLimitMin(c.fBckFitLimitMin),
123 fBckFitLimitMax(c.fBckFitLimitMax),
124 fInvariantMassFit(c.fInvariantMassFit),
125 fkAnalysisType(c.fkAnalysisType),
133 //___________________________________________________________________________
134 AliAnalysisTaskMuonDistributions::~AliAnalysisTaskMuonDistributions() {
138 Info("~AliAnalysisTaskMuonDistributions","Calling Destructor");
146 //___________________________________________________________________________
147 void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
149 // output objects creation
151 fOutput = new TList();
156 TH1D *hNumberMuonTracks = new TH1D("hNumberMuonTracks","hNumberMuonTracks;N_{#mu tracks}",10,0.,10.);
160 TH1D *hMassDimu = new TH1D("hMassDimu","hMassDimu;M_{#mu#mu} (GeV/c^{2})",180,0,9.);
161 TH1D *hPtDimu = new TH1D("hPtDimu","hPtDimu;p_{T} (GeV/c)",100,0,20);
162 TH1D *hRapidityDimu = new TH1D("hRapidityDimu","hRapidityDimu;y",100,-5,-2);
163 TH1D *hCosThetaCSDimu = new TH1D("hCosThetaCSDimu","hCosThetaCSDimu;cos#theta_{CS}",100,-1.,1.);
164 TH1D *hCosThetaHEDimu = new TH1D("hCosThetaHEDimu","hCosThetaHEDimu;cos#theta_{HE}",100,-1.,1.);
168 TH1D *hP = new TH1D("hP","hP;p (GeV/c)",100,0,500);
169 TH1D *hPt = new TH1D("hPt","hPt;p_{T} (GeV/c)",100,0,20);
170 TH1D *hRapidity = new TH1D("hRapidity","hRapidity;y",100,-5,-2);
172 fOutput->Add(hNumberMuonTracks);
173 fOutput->Add(hMassDimu);
174 fOutput->Add(hPtDimu);
175 fOutput->Add(hRapidityDimu);
176 fOutput->Add(hCosThetaCSDimu);
177 fOutput->Add(hCosThetaHEDimu);
180 fOutput->Add(hRapidity);
184 //_________________________________________________
185 void AliAnalysisTaskMuonDistributions::UserExec(Option_t *)
188 // Execute analysis for current event
190 AliESDEvent *esd=0x0;
191 AliAODEvent *aod=0x0;
193 if(strcmp(fkAnalysisType,"ESD")==0){
194 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
195 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
196 esd = esdH->GetEvent();
197 } else if(strcmp(fkAnalysisType,"AOD")==0){
198 aod = dynamic_cast<AliAODEvent*> (InputEvent());
202 if(strcmp(fkAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks();
203 else if(strcmp(fkAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks();
206 for (Int_t j = 0; j<ntracks; j++) {
207 Float_t pxmu1=-999; Float_t pymu1=-999; Float_t pzmu1=-999; Float_t ptmu1=-999; Float_t pmu1=-999;
208 Float_t emu1=-999; Float_t rapiditymu1=-999; Float_t chargemu1=-999;
209 if(strcmp(fkAnalysisType,"ESD")==0){
210 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(esd->GetMuonTrack(j)));
211 if (!mu1->ContainTrackerData()) continue;
212 chargemu1 = mu1->Charge();
219 rapiditymu1 = Rapidity(emu1,pzmu1);
220 } else if(strcmp(fkAnalysisType,"AOD")==0){
221 AliAODTrack *mu1 = aod->GetTrack(j);
222 if(!mu1->IsMuonTrack()) continue;
223 chargemu1 = mu1->Charge();
230 rapiditymu1 = Rapidity(emu1,pzmu1);
232 ((TH1D*)(fOutput->FindObject("hP")))->Fill(pmu1);
233 ((TH1D*)(fOutput->FindObject("hPt")))->Fill(ptmu1);
234 ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapiditymu1);
237 for (Int_t jj = 0; jj<ntracks; jj++) {
238 Float_t pxmu2=-999; Float_t pymu2=-999; Float_t pzmu2=-999;
239 Float_t emu2=-999;Float_t chargemu2=-999;
240 if(strcmp(fkAnalysisType,"ESD")==0){
241 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(esd->GetMuonTrack(jj)));
242 if (!mu2->ContainTrackerData()) continue;
243 chargemu2 = mu2->Charge();
248 } else if(strcmp(fkAnalysisType,"AOD")==0){
249 AliAODTrack *mu2 = aod->GetTrack(jj);
250 if(!mu2->IsMuonTrack()) continue;
251 chargemu2 = mu2->Charge();
258 Float_t ptdimu = TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
259 Float_t massdimu = InvMass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
260 Float_t rapiditydimu = Rapidity((emu1+emu2),(pzmu1+pzmu2));
261 Double_t costhetaCSdimu = CostCS((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
262 Double_t costhetaHEdimu = CostHE((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
263 ((TH1D*)(fOutput->FindObject("hMassDimu")))->Fill(massdimu);
264 ((TH1D*)(fOutput->FindObject("hPtDimu")))->Fill(ptdimu);
265 ((TH1D*)(fOutput->FindObject("hRapidityDimu")))->Fill(rapiditydimu);
266 ((TH1D*)(fOutput->FindObject("hCosThetaCSDimu")))->Fill(costhetaCSdimu);
267 ((TH1D*)(fOutput->FindObject("hCosThetaHEDimu")))->Fill(costhetaHEdimu);
274 ((TH1D*)(fOutput->FindObject("hNumberMuonTracks")))->Fill(nmuontracks);
280 //________________________________________________________________________
281 void AliAnalysisTaskMuonDistributions::Terminate(Option_t *)
286 gStyle->SetCanvasColor(10);
287 gStyle->SetFrameFillColor(10);
291 printf("Using beam Energy=%f \n",fBeamEnergy);
293 fOutput = dynamic_cast<TList*> (GetOutputData(1));
294 TH1D *hNumberMuonTracks = dynamic_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));
295 TH1D *hMassDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hMassDimu"));
296 TH1D *hPtDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hPtDimu"));
297 TH1D *hRapidityDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidityDimu"));
298 TH1D *hCosThetaCSDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));
299 TH1D *hCosThetaHEDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));
300 TH1D *hP = dynamic_cast<TH1D*> (fOutput->FindObject("hP"));
301 TH1D *hPt = dynamic_cast<TH1D*> (fOutput->FindObject("hPt"));
302 TH1D *hRapidity = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidity"));
304 TCanvas *c0_MuonDistributions = new TCanvas("c0_MuonDistributions","Plots",xmin,ymin,600,600);
305 c0_MuonDistributions->Divide(2,2);
306 c0_MuonDistributions->cd(1);
307 hNumberMuonTracks->SetFillColor(2);
308 hNumberMuonTracks->Draw();
311 TCanvas *c1_MuonDistributions = new TCanvas("c1_MuonDistributions","Muon kinematic distributions Plots",xmin,ymin,600,600);
312 c1_MuonDistributions->Divide(2,2);
313 c1_MuonDistributions->cd(1);
317 c1_MuonDistributions->cd(2);
319 hPt->SetLineColor(4);
321 c1_MuonDistributions->cd(3);
322 hRapidity->SetLineColor(4);
326 TCanvas *c2_MuonDistributions = new TCanvas("c2_MuonDistributions","Dimuon kinematic distributions Plots",xmin,ymin,600,600);
327 c2_MuonDistributions->Divide(2,2);
328 c2_MuonDistributions->cd(1);
330 hPtDimu->SetLineColor(2);
332 c2_MuonDistributions->cd(2);
333 hRapidityDimu->SetLineColor(2);
334 hRapidityDimu->Draw();
335 c2_MuonDistributions->cd(3);
336 hCosThetaCSDimu->SetLineColor(2);
337 hCosThetaCSDimu->Draw();
338 c2_MuonDistributions->cd(4);
339 hCosThetaHEDimu->SetLineColor(2);
340 hCosThetaHEDimu->Draw();
343 TCanvas *c3_MuonDistributions = new TCanvas("c3_MuonDistributions","Invariant Mass Plots",xmin,ymin,600,600);
346 if(fInvariantMassFit && hMassDimu->GetEntries()>100) FitInvMass(hMassDimu);
347 c3_MuonDistributions->Update();
350 //________________________________________________________________________
351 Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
352 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
355 // invariant mass calculation
357 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
358 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
361 //________________________________________________________________________
362 Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const
365 // calculate rapidity
368 if(TMath::Abs(e-pz)>1e-7){
369 rap = 0.5*TMath::Log((e+pz)/(e-pz));
376 //________________________________________________________________________
377 Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
378 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
381 // costCS calculation
383 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
384 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
385 TVector3 beta,zaxisCS;
386 Double_t mp=0.93827231;
388 // --- Fill the Lorentz vector for projectile and target in the CM frame
390 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
391 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
393 // --- Get the muons parameters in the CM frame
395 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
396 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
398 // --- Obtain the dimuon parameters in the CM frame
400 pDimuCM=pMu1CM+pMu2CM;
402 // --- Translate the dimuon parameters in the dimuon rest frame
404 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
409 pMu1Dimu.Boost(beta);
410 pMu2Dimu.Boost(beta);
411 pProjDimu.Boost(beta);
412 pTargDimu.Boost(beta);
414 // --- Determine the z axis for the CS angle
416 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
418 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
421 if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
422 else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
426 //________________________________________________________________________
427 Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
428 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
431 // costHE calculation
433 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
434 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
435 TVector3 beta,zaxisCS;
437 // --- Get the muons parameters in the CM frame
439 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
440 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
442 // --- Obtain the dimuon parameters in the CM frame
444 pDimuCM=pMu1CM+pMu2CM;
446 // --- Translate the muon parameters in the dimuon rest frame
448 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
451 pMu1Dimu.Boost(beta);
452 pMu2Dimu.Boost(beta);
454 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
457 zaxis=(pDimuCM.Vect()).Unit();
459 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
461 if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
462 else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
466 //________________________________________________________________________
467 void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
470 // Fit to the Invariant Mass Spectrum
472 TF1 *gaupsi = new TF1("gaupsi","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
473 TF1 *gaupsip = new TF1("gaupsip","gaus",fPsiPFitLimitMin,fPsiPFitLimitMax);
474 TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);
475 TF1 *tot = new TF1("mtot","gaus(0)+expo(3)+gaus(5)",fInvMassFitLimitMin,fInvMassFitLimitMax);
477 Double_t binWidth= histo->GetBinWidth(1);
478 gaupsi->SetLineColor(kGreen);
479 gaupsi->SetLineWidth(2);
480 histo->Fit(gaupsi,"Rl0");
481 ex->SetLineColor(kBlue);
483 histo->Fit(ex,"Rl+");
484 gaupsip->SetLineColor(kAzure-9);
485 gaupsip->SetLineWidth(2);
486 gaupsip->SetParLimits(1,3.6,3.75);
487 gaupsip->SetParLimits(2,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2));
488 histo->Fit(gaupsip,"Rl0+");
489 gaupsi->GetParameters(&par[0]);
490 ex->GetParameters(&par[3]);
491 gaupsip->GetParameters(&par[5]);
492 tot->SetParameters(par);
493 tot->SetLineColor(2);
494 tot->SetLineWidth(2);
495 tot->SetParLimits(6,3.6,3.75); //limit for the psi(2S) parameters
496 tot->SetParLimits(7,0.8*gaupsi->GetParameter(2),1.2*gaupsi->GetParameter(2)); //limit for the psi(2S) parameters
497 histo->Fit(tot,"Rl+");
499 Double_t chi2 = tot->GetChisquare();
500 Double_t ndf = tot->GetNDF();
502 Float_t meanPsi= tot->GetParameter(1);
503 Float_t sigPsi= tot->GetParameter(2)*1000.;
504 Double_t nPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/binWidth;
505 TF1 *psifix = new TF1("psifix","gaus",2.,5.);
506 psifix->SetParameter(0,tot->GetParameter(0));
507 psifix->SetParameter(1,tot->GetParameter(1));
508 psifix->SetParameter(2,tot->GetParameter(2));
509 psifix->SetLineColor(kGreen);
510 psifix->Draw("same");
511 Double_t nPsi2933 = psifix->Integral(2.9,3.3)/binWidth;
513 TF1 *exfix = new TF1("exfix","expo",2.,5.);
514 exfix->SetParameter(0,tot->GetParameter(3));
515 exfix->SetParameter(1,tot->GetParameter(4));
516 Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
518 Float_t meanPsiP= tot->GetParameter(6);
519 Float_t sigPsiP= tot->GetParameter(7)*1000.;
520 Double_t nPsiPFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(5)*tot->GetParameter(7)/binWidth;
521 TF1 *psipfix = new TF1("psipfix","gaus",3.,5.);
522 psipfix->SetParameter(0,tot->GetParameter(5));
523 psipfix->SetParameter(1,tot->GetParameter(6));
524 psipfix->SetParameter(2,tot->GetParameter(7));
525 psipfix->SetLineColor(kAzure-9);
526 psipfix->Draw("same");
528 printf("\n\n****************************************************************************\n");
530 if(nPsiFit<0) nPsiFit = 0;
531 sprintf(psitext,"N. J/#psi = %10.0f",nPsiFit);
532 printf("\nN. J/psi = %10.0f\n",nPsiFit);
533 TLatex *psilatex = new TLatex(4.,0.85*histo->GetMaximum(),psitext);
534 psilatex->SetTextColor(2);
535 psilatex->SetTextSize(0.03);
536 psilatex->SetTextAlign(2);
540 sprintf(psi2text,"J/#psi m=%4.3f GeV #sigma= %4.2f MeV",meanPsi,sigPsi);
541 printf("J/psi m= %4.3f GeV sigma= %4.2f MeV\n",meanPsi,sigPsi);
542 TLatex *psi2latex = new TLatex(4.,0.425*histo->GetMaximum(),psi2text);
543 psi2latex->SetTextColor(2);
544 psi2latex->SetTextSize(0.03);
545 psi2latex->SetTextAlign(2);
549 sprintf(sbtext,"S/B (2.9-3.3)= %4.2f ",nPsi2933/nBck);
550 printf("S/B (2.9-3.3)= %4.2f\n",nPsi2933/nBck);
551 TLatex *sblatex = new TLatex(4.,0.212*histo->GetMaximum(),sbtext);
552 sblatex->SetTextColor(2);
553 sblatex->SetTextSize(0.03);
554 sblatex->SetTextAlign(2);
558 if(nPsiPFit<0) nPsiPFit = 0;
559 sprintf(psiptext,"N. #psi(2S) = %10.0f",nPsiPFit);
560 printf("\npsi(2S) = %10.0f\n",nPsiPFit);
561 TLatex *psiplatex = new TLatex(4.,0.106*histo->GetMaximum(),psiptext);
562 psiplatex->SetTextColor(2);
563 psiplatex->SetTextSize(0.03);
564 psiplatex->SetTextAlign(2);
568 sprintf(psip2text,"#psi(2S) m=%4.3f GeV #sigma= %4.2f MeV",meanPsiP,sigPsiP);
569 printf("psi(2S) m= %4.3f GeV sigma= %4.2f MeV\n",meanPsiP,sigPsiP);
570 TLatex *psip2latex = new TLatex(4.,0.053*histo->GetMaximum(),psip2text);
571 psip2latex->SetTextColor(2);
572 psip2latex->SetTextSize(0.03);
573 psip2latex->SetTextAlign(2);
577 sprintf(chi2text,"#chi^2/ndf = %4.2f ",chi2/ndf);
578 printf("chi^2/ndf = %4.2f\n",chi2/ndf);
579 TLatex *chi2latex = new TLatex(4.,0.026*histo->GetMaximum(),chi2text);
580 chi2latex->SetTextColor(2);
581 chi2latex->SetTextSize(0.03);
582 chi2latex->SetTextAlign(2);
584 printf("\n****************************************************************************\n");