From e4dddba1aab57920bc1d55c8bc94fda885da7410 Mon Sep 17 00:00:00 2001 From: martinez Date: Fri, 9 Oct 2009 07:25:22 +0000 Subject: [PATCH] Analysis of muon and dimuon kinematical distribution for the official analysis train ESD->AOD (Roberta) --- PWG3/CMake_libPWG3muon.txt | 1 + PWG3/PWG3muonLinkDef.h | 1 + PWG3/libPWG3muon.pkg | 3 +- PWG3/muon/AddTaskMuonDistributions.C | 44 ++ .../muon/AliAnalysisTaskMuonDistributions.cxx | 507 ++++++++++++++++++ PWG3/muon/AliAnalysisTaskMuonDistributions.h | 59 ++ 6 files changed, 614 insertions(+), 1 deletion(-) create mode 100644 PWG3/muon/AddTaskMuonDistributions.C create mode 100755 PWG3/muon/AliAnalysisTaskMuonDistributions.cxx create mode 100755 PWG3/muon/AliAnalysisTaskMuonDistributions.h diff --git a/PWG3/CMake_libPWG3muon.txt b/PWG3/CMake_libPWG3muon.txt index 5ee8273f63e..2bebf99e3bd 100644 --- a/PWG3/CMake_libPWG3muon.txt +++ b/PWG3/CMake_libPWG3muon.txt @@ -16,6 +16,7 @@ set(SRCS muon/AliEventPoolMuon.cxx muon/AliAnalysisTaskCreateMixedDimuons.cxx muon/AliAnalysisTaskMuonAODCreation.cxx + muon/AliAnalysisTaskMuonDistributions.cxx ) # fill list of header files from list of source files diff --git a/PWG3/PWG3muonLinkDef.h b/PWG3/PWG3muonLinkDef.h index c440df88cbe..82eefb34914 100644 --- a/PWG3/PWG3muonLinkDef.h +++ b/PWG3/PWG3muonLinkDef.h @@ -19,6 +19,7 @@ #pragma link C++ class AliEventPoolMuon+; #pragma link C++ class AliAnalysisTaskCreateMixedDimuons+; #pragma link C++ class AliAnalysisTaskMuonAODCreation+; +#pragma link C++ class AliAnalysisTaskMuonDistributions+; #endif diff --git a/PWG3/libPWG3muon.pkg b/PWG3/libPWG3muon.pkg index f5417c30572..fbb5a8ffe59 100644 --- a/PWG3/libPWG3muon.pkg +++ b/PWG3/libPWG3muon.pkg @@ -14,7 +14,8 @@ SRCS:= muon/AliAnalysisTaskESDMuonFilter.cxx \ muon/AliCFMuonResTask1.cxx \ muon/AliEventPoolMuon.cxx \ muon/AliAnalysisTaskCreateMixedDimuons.cxx \ - muon/AliAnalysisTaskMuonAODCreation.cxx + muon/AliAnalysisTaskMuonAODCreation.cxx \ + muon/AliAnalysisTaskMuonDistributions.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/PWG3/muon/AddTaskMuonDistributions.C b/PWG3/muon/AddTaskMuonDistributions.C new file mode 100644 index 00000000000..b90fe75a8f5 --- /dev/null +++ b/PWG3/muon/AddTaskMuonDistributions.C @@ -0,0 +1,44 @@ +AliAnalysisTaskMuonDistributions *AddTaskMuonDistributions(const char *kAnalysisType){ + +//**************************************************************************************** +// Add task class. +// The attached class prepares and draws some kinematical distributions of muons/dimuons +// Roberta +//**************************************************************************************** + + printf("Creating Task for Muon/Dimuon Histos\n"); + + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskMuonDistributions", "No analysis manager to connect to."); + return NULL; + } + + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist1",TList::Class(),AliAnalysisManager::kOutputContainer,"MuonDistributions.root"); + + AliAnalysisTaskMuonDistributions *MuonDistributionsTask = new AliAnalysisTaskMuonDistributions("AliAnalysisTaskMuonDistributions"); + MuonDistributionsTask->SetAnalysisType(kAnalysisType); + // + // define by hand the beam energy + // + MuonDistributionsTask->SetBeamEnergy(5000.); + // + // define fits limits + // + MuonDistributionsTask->SetInvMassFitLimits(2.,5.5); + MuonDistributionsTask->SetPsiFitLimits(2.9,3.3); + MuonDistributionsTask->SetBckFitLimits(2.,2.8); + // + // perform fit to the invariant mass spectrum + // + MuonDistributionsTask->FitInvariantMassSpectrum(kFALSE); + + mgr->AddTask(MuonDistributionsTask); + + mgr->ConnectInput(MuonDistributionsTask,0,mgr->GetCommonInputContainer()); + mgr->ConnectOutput(MuonDistributionsTask,1,coutput1); + + return MuonDistributionsTask; +} diff --git a/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx b/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx new file mode 100755 index 00000000000..6fc32de980a --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx @@ -0,0 +1,507 @@ +/************************************************************************** + * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//----------------------------------------------------------------------------- +// Analysis task to compute muon/dimuon kinematic distributions +// The output is a list of histograms. +// The macro class can run on AOD or in the train with the ESD filter. +// R. Arnaldi +// +//----------------------------------------------------------------------------- + +//#ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_CXX +//#define ALIANALYSISTASKMUONDISTRIBUTIONS_CXX + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliAnalysisTaskMuonDistributions.h" +#include "AliAnalysisTaskSE.h" +#include "AliAnalysisManager.h" +#include "AliAnalysisDataSlot.h" +#include "AliESDEvent.h" +#include "AliESD.h" +#include "AliVEvent.h" +#include "AliMCEventHandler.h" +#include "AliInputEventHandler.h" +#include "AliMCEvent.h" +#include "AliStack.h" +#include "AliLog.h" +#include "AliHeader.h" +#include "AliESDHeader.h" +#include "AliStack.h" +#include "TParticle.h" +#include "TLorentzVector.h" +#include "AliESDMuonTrack.h" +#include "AliESDtrack.h" +#include "AliESDInputHandler.h" +#include "AliAODEvent.h" +#include "AliAODHeader.h" +#include "AliAODHandler.h" +#include "AliAODInputHandler.h" +#include "AliMCEventHandler.h" + +ClassImp(AliAnalysisTaskMuonDistributions) + +//__________________________________________________________________________ +AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions() : + fBeamEnergy(0.), + fInvMassFitLimitMin(2.), + fInvMassFitLimitMax(5.), + fPsiFitLimitMin(2.9), + fPsiFitLimitMax(3.3), + fBckFitLimitMin(2.2), + fBckFitLimitMax(2.85), + fInvariantMassFit(kFALSE), + fAnalysisType(0x0), + fOutput(0x0) +{ +} +//___________________________________________________________________________ +AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const Char_t* name) : + AliAnalysisTaskSE(name), + fBeamEnergy(0.), + fInvMassFitLimitMin(2.), + fInvMassFitLimitMax(5.), + fPsiFitLimitMin(2.9), + fPsiFitLimitMax(3.3), + fBckFitLimitMin(2.2), + fBckFitLimitMax(2.85), + fInvariantMassFit(kFALSE), + fAnalysisType(0x0), + fOutput(0x0) +{ + // Constructor. Initialization of Inputs and Outputs + // + Info("AliAnalysisTaskMuonDistributions","Calling Constructor"); + + DefineOutput(1,TList::Class()); + +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonDistributions& AliAnalysisTaskMuonDistributions::operator=(const AliAnalysisTaskMuonDistributions& c) +{ + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c) ; + } + return *this; +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const AliAnalysisTaskMuonDistributions& c) : + AliAnalysisTaskSE(c), + fBeamEnergy(c.fBeamEnergy), + fInvMassFitLimitMin(c.fInvMassFitLimitMin), + fInvMassFitLimitMax(c.fInvMassFitLimitMax), + fPsiFitLimitMin(c.fPsiFitLimitMin), + fPsiFitLimitMax(c.fPsiFitLimitMax), + fBckFitLimitMin(c.fBckFitLimitMin), + fBckFitLimitMax(c.fBckFitLimitMax), + fInvariantMassFit(c.fInvariantMassFit), + fAnalysisType(c.fAnalysisType), + fOutput(c.fOutput) + { + // + // Copy Constructor + // +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonDistributions::~AliAnalysisTaskMuonDistributions() { + // + //destructor + // + Info("~AliAnalysisTaskMuonDistributions","Calling Destructor"); +} + +//___________________________________________________________________________ +void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){ + + fOutput = new TList(); + fOutput->SetOwner(); + // + // various histos + // + TH1D *hNumberMuonTracks = new TH1D("hNumberMuonTracks","hNumberMuonTracks;N_{#mu tracks}",10,0.,10.); + // + // dimuon histos + // + TH1D *hMass_Dimu = new TH1D("hMass_Dimu","hMass_Dimu;M_{#mu#mu} (GeV/c^{2})",180,0,9.); + TH1D *hPt_Dimu = new TH1D("hPt_Dimu","hPt_Dimu;p_{T} (GeV/c)",100,0,20); + TH1D *hRapidity_Dimu = new TH1D("hRapidity_Dimu","hRapidity_Dimu;y",100,-5,-2); + TH1D *hCosThetaCS_Dimu = new TH1D("hCosThetaCS_Dimu","hCosThetaCS_Dimu;cos#theta_{CS}",100,-1.,1.); + TH1D *hCosThetaHE_Dimu = new TH1D("hCosThetaHE_Dimu","hCosThetaHE_Dimu;cos#theta_{HE}",100,-1.,1.); + // + // muon histos + // + TH1D *hP = new TH1D("hP","hP;p (GeV/c)",100,0,500); + TH1D *hPt = new TH1D("hPt","hPt;p_{T} (GeV/c)",100,0,20); + TH1D *hRapidity = new TH1D("hRapidity","hRapidity;y",100,-5,-2); + + fOutput->Add(hNumberMuonTracks); + fOutput->Add(hMass_Dimu); + fOutput->Add(hPt_Dimu); + fOutput->Add(hRapidity_Dimu); + fOutput->Add(hCosThetaCS_Dimu); + fOutput->Add(hCosThetaHE_Dimu); + fOutput->Add(hP); + fOutput->Add(hPt); + fOutput->Add(hRapidity); + fOutput->ls(); +} + +//_________________________________________________ +void AliAnalysisTaskMuonDistributions::UserExec(Option_t *) +{ + AliESDEvent *esd=0x0; + AliAODEvent *aod=0x0; + + if(strcmp(fAnalysisType,"ESD")==0){ + AliESDInputHandler *esdH = dynamic_cast + (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + esd = esdH->GetEvent(); + } else if(strcmp(fAnalysisType,"AOD")==0){ + aod = dynamic_cast (InputEvent()); + } + + Int_t ntracks=-999; + if(strcmp(fAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks(); + else if(strcmp(fAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks(); + Int_t nmuontracks=0; + + for (Int_t j = 0; jGetMuonTrack(j))); + if (!mu1->ContainTrackerData()) continue; + charge_mu1 = mu1->Charge(); + px_mu1 = mu1->Px(); + py_mu1 = mu1->Py(); + pz_mu1 = mu1->Pz(); + e_mu1 = mu1->E(); + p_mu1 = mu1->P(); + pt_mu1 = mu1->Pt(); + rapidity_mu1 = Rapidity(e_mu1,pz_mu1); + } else if(strcmp(fAnalysisType,"AOD")==0){ + AliAODTrack *mu1 = aod->GetTrack(j); + if(!mu1->IsMuonTrack()) continue; + charge_mu1 = mu1->Charge(); + px_mu1 = mu1->Px(); + py_mu1 = mu1->Py(); + pz_mu1 = mu1->Pz(); + e_mu1 = mu1->E(); + p_mu1 = mu1->P(); + pt_mu1 = mu1->Pt(); + rapidity_mu1 = Rapidity(e_mu1,pz_mu1); + } + ((TH1D*)(fOutput->FindObject("hP")))->Fill(p_mu1); + ((TH1D*)(fOutput->FindObject("hPt")))->Fill(pt_mu1); + ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapidity_mu1); + nmuontracks++; + if(charge_mu1<0){ + for (Int_t jj = 0; jjGetMuonTrack(jj))); + if (!mu2->ContainTrackerData()) continue; + charge_mu2 = mu2->Charge(); + px_mu2 = mu2->Px(); + py_mu2 = mu2->Py(); + pz_mu2 = mu2->Pz(); + e_mu2 = mu2->E(); + } else if(strcmp(fAnalysisType,"AOD")==0){ + AliAODTrack *mu2 = aod->GetTrack(jj); + if(!mu2->IsMuonTrack()) continue; + charge_mu2 = mu2->Charge(); + px_mu2 = mu2->Px(); + py_mu2 = mu2->Py(); + pz_mu2 = mu2->Pz(); + e_mu2 = mu2->E(); + } + if(charge_mu2>0){ + Float_t pt_dimu = TMath::Sqrt((px_mu1+px_mu2)*(px_mu1+px_mu2)+(py_mu1+py_mu2)*(py_mu1+py_mu2)); + Float_t mass_dimu = InvMass(e_mu1,px_mu1,py_mu1,pz_mu1,e_mu2,px_mu2,py_mu2,pz_mu2); + Float_t rapidity_dimu = Rapidity((e_mu1+e_mu2),(pz_mu1+pz_mu2)); + Double_t costhetaCS_dimu = CostCS((Double_t) px_mu1,(Double_t) py_mu1,(Double_t)pz_mu1,(Double_t) e_mu1,(Double_t)charge_mu1,(Double_t) px_mu2,(Double_t) py_mu2,(Double_t)pz_mu2,(Double_t) e_mu2); + Double_t costhetaHE_dimu = CostHE((Double_t) px_mu1,(Double_t) py_mu1,(Double_t)pz_mu1,(Double_t) e_mu1,(Double_t)charge_mu1,(Double_t) px_mu2,(Double_t) py_mu2,(Double_t)pz_mu2,(Double_t) e_mu2); + ((TH1D*)(fOutput->FindObject("hMass_Dimu")))->Fill(mass_dimu); + ((TH1D*)(fOutput->FindObject("hPt_Dimu")))->Fill(pt_dimu); + ((TH1D*)(fOutput->FindObject("hRapidity_Dimu")))->Fill(rapidity_dimu); + ((TH1D*)(fOutput->FindObject("hCosThetaCS_Dimu")))->Fill(costhetaCS_dimu); + ((TH1D*)(fOutput->FindObject("hCosThetaHE_Dimu")))->Fill(costhetaHE_dimu); + } + //delete mu2; + } // second mu Loop + } // mu- Selection + //delete mu1; + } + ((TH1D*)(fOutput->FindObject("hNumberMuonTracks")))->Fill(nmuontracks); + + PostData(1,fOutput); + } + + +//________________________________________________________________________ +void AliAnalysisTaskMuonDistributions::Terminate(Option_t *) +{ + gStyle->SetCanvasColor(10); + gStyle->SetFrameFillColor(10); + Int_t xmin=20; + Int_t ymin=20; + + printf("Using beam Energy=%f \n",fBeamEnergy); + + TH1D *h_NumberMuonTracks = dynamic_cast (fOutput->FindObject("hNumberMuonTracks")); + TH1D *h_Mass_Dimu = dynamic_cast (fOutput->FindObject("hMass_Dimu")); + TH1D *h_Pt_Dimu = dynamic_cast (fOutput->FindObject("hPt_Dimu")); + TH1D *h_Rapidity_Dimu = dynamic_cast (fOutput->FindObject("hRapidity_Dimu")); + TH1D *h_CostCS_Dimu = dynamic_cast (fOutput->FindObject("hCosThetaCS_Dimu")); + TH1D *h_CostHE_Dimu = dynamic_cast (fOutput->FindObject("hCosThetaHE_Dimu")); + TH1D *h_P = dynamic_cast (fOutput->FindObject("hP")); + TH1D *h_Pt = dynamic_cast (fOutput->FindObject("hPt")); + TH1D *h_Rapidity = dynamic_cast (fOutput->FindObject("hRapidity")); + + TCanvas *c0 = new TCanvas("c0","Plots",xmin,ymin,600,600); + c0->Divide(2,2); + c0->cd(1); + h_NumberMuonTracks->Draw(); + + xmin+=20; ymin+=20; + TCanvas *c1 = new TCanvas("c1","Muon kinematic distributions Plots",xmin,ymin,600,600); + c1->Divide(2,2); + c1->cd(1); + gPad->SetLogy(1); + h_P->Draw(); + c1->cd(2); + gPad->SetLogy(1); + h_Pt->Draw(); + c1->cd(3); + h_Rapidity->Draw(); + + xmin+=20; ymin+=20; + TCanvas *c2 = new TCanvas("c2","Dimuon kinematic distributions Plots",xmin,ymin,600,600); + c2->Divide(2,2); + c2->cd(1); + gPad->SetLogy(1); + h_Pt_Dimu->Draw(); + c2->cd(2); + h_Rapidity_Dimu->Draw(); + c2->cd(3); + h_CostCS_Dimu->Draw(); + c2->cd(4); + h_CostHE_Dimu->Draw(); + + xmin+=20; ymin+=20; + TCanvas *c3 = new TCanvas("c3","Invariant Mass Plots",xmin,ymin,600,600); + gPad->SetLogy(1); + h_Mass_Dimu->Draw(); + if(fInvariantMassFit) FitInvMass(h_Mass_Dimu); + c3->Update(); + +} + +//________________________________________________________________________ +Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1, + Float_t e2, Float_t px2, Float_t py2, Float_t pz2) +{ +// invariant mass calculation + Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+ + (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2))); + return imassrec; +} +//________________________________________________________________________ +Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) +{ +// calculate rapidity + Float_t rap; + if(e!=pz){ + rap = 0.5*TMath::Log((e+pz)/(e-pz)); + return rap; + } else { + rap = -200; + return rap; + } +} +//________________________________________________________________________ +Double_t AliAnalysisTaskMuonDistributions::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2) +{ + TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,zaxisCS; + Double_t mp=0.93827231; + // + // --- Fill the Lorentz vector for projectile and target in the CM frame + // + PProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); + PTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the dimuon parameters in the dimuon rest frame + // + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PProjDimu=PProjCM; + PTargDimu=PTargCM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + PProjDimu.Boost(beta); + PTargDimu.Boost(beta); + // + // --- Determine the z axis for the CS angle + // + zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit(); + // + // --- Determine the CS angle (angle between mu+ and the z axis defined above) + Double_t cost; + + if(charge1>0) {cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());} + else {cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());} + return cost; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskMuonDistributions::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2) +{ + TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame + TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame + TVector3 beta,zaxisCS; + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the muon parameters in the dimuon rest frame + // + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + // + // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system) + // + TVector3 zaxis; + zaxis=(PDimuCM.Vect()).Unit(); + + // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above) + Double_t cost; + if(charge1>0) {cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());} + else {cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());} + return cost; +} + +//________________________________________________________________________ +void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo) +{ + TF1 *gau = new TF1("gau","gaus",fPsiFitLimitMin,fPsiFitLimitMax); + TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax); + TF1 *tot = new TF1("mtot","gaus(0)+expo(3)",fInvMassFitLimitMin,fInvMassFitLimitMax); + Double_t par[5]; + Double_t BinWidth= histo->GetBinWidth(1); + printf("BW=%f\n",BinWidth); + gau->SetLineColor(3); + gau->SetLineWidth(2); + histo->Fit(gau,"RlQ"); + ex->SetLineColor(4); + ex->SetLineWidth(2); + histo->Fit(ex,"RlQ+"); + gau->GetParameters(&par[0]); + ex->GetParameters(&par[3]); + tot->SetParameters(par); + tot->SetLineColor(2); + tot->SetLineWidth(2); + histo->Fit(tot,"Rl+"); + histo->Draw("e"); + Double_t Chi2 = tot->GetChisquare(); + Double_t NDF = tot->GetNDF(); + Float_t MeanPsi= tot->GetParameter(1); + Float_t SigPsi= tot->GetParameter(2)*1000.; + Double_t NPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/BinWidth; + TF1 *exfix = new TF1("exfix","expo",2.,5.); + exfix->SetParameter(0,tot->GetParameter(3)); + exfix->SetParameter(1,tot->GetParameter(4)); + Double_t NBck = exfix->Integral(2.9,3.3)/BinWidth; + + printf("\n\n****************************************************************************\n"); + char psi_text[100]; + sprintf(psi_text,"N. J/#psi (2.9-3.3)=%10.0f",NPsiFit); + printf("\nN. J/psi (2.9-3.3)=%10.0f\n",NPsiFit); + TLatex *psi_latex = new TLatex(4.5,0.85*histo->GetMaximum(),psi_text); + psi_latex->SetTextColor(2); + psi_latex->SetTextSize(0.03); + psi_latex->SetTextAlign(2); + psi_latex->Draw(); + + char psi2_text[100]; + sprintf(psi2_text,"J/#psi m=%4.3f GeV #sigma=%4.2f MeV",MeanPsi,SigPsi); + printf("J/psi m=%4.3f GeV sigma=%4.2f MeV\n",MeanPsi,SigPsi); + TLatex *psi2_latex = new TLatex(4.5,0.425*histo->GetMaximum(),psi2_text); + psi2_latex->SetTextColor(2); + psi2_latex->SetTextSize(0.03); + psi2_latex->SetTextAlign(2); + psi2_latex->Draw(); + + char sb_text[100]; + sprintf(sb_text,"S/B (2.9-3.3)=%4.2f ",NPsiFit/NBck); + printf("S/B (2.9-3.3)=%4.2f\n",NPsiFit/NBck); + TLatex *sb_latex = new TLatex(4.5,0.212*histo->GetMaximum(),sb_text); + sb_latex->SetTextColor(2); + sb_latex->SetTextSize(0.03); + sb_latex->SetTextAlign(2); + sb_latex->Draw(); + + char chi2_text[100]; + sprintf(chi2_text,"#chi^2/ndf =%4.2f ",Chi2/NDF); + printf("chi^2/ndf =%4.2f\n",Chi2/NDF); + TLatex *chi2_latex = new TLatex(4.5,0.106*histo->GetMaximum(),chi2_text); + chi2_latex->SetTextColor(2); + chi2_latex->SetTextSize(0.03); + chi2_latex->SetTextAlign(2); + chi2_latex->Draw(); + printf("\n****************************************************************************\n"); + +} + + diff --git a/PWG3/muon/AliAnalysisTaskMuonDistributions.h b/PWG3/muon/AliAnalysisTaskMuonDistributions.h new file mode 100755 index 00000000000..c485e12b754 --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonDistributions.h @@ -0,0 +1,59 @@ +#ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_H +#define ALIANALYSISTASKMUONDISTRIBUTIONS_H + +#include "AliAnalysisTaskSE.h" +#include "TMath.h" + +class TH1D; +class TParticle ; +class TLorentzVector ; +class TFile ; +class AliStack ; +class AliESDtrack; +class AliVParticle; + +class AliAnalysisTaskMuonDistributions : public AliAnalysisTaskSE { + public: + + AliAnalysisTaskMuonDistributions(); + AliAnalysisTaskMuonDistributions(const Char_t* name); + AliAnalysisTaskMuonDistributions& operator= (const AliAnalysisTaskMuonDistributions& c); + AliAnalysisTaskMuonDistributions(const AliAnalysisTaskMuonDistributions& c); + virtual ~AliAnalysisTaskMuonDistributions(); + + // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects + void UserExec(Option_t *option); + void Terminate(Option_t *); + void UserCreateOutputObjects(); + + void SetBeamEnergy(Double_t en) {fBeamEnergy=en;} + void SetAnalysisType(const char* type) {fAnalysisType=type;} + void SetInvMassFitLimits(Double_t xmin, Double_t xmax) {fInvMassFitLimitMin=xmin; fInvMassFitLimitMax=xmax;} + void SetPsiFitLimits(Double_t xmin, Double_t xmax) {fPsiFitLimitMin=xmin; fPsiFitLimitMax=xmax;} + void SetBckFitLimits(Double_t xmin, Double_t xmax) {fBckFitLimitMin=xmin; fBckFitLimitMax=xmax;} + void FitInvariantMassSpectrum(Bool_t massfit=kFALSE) {fInvariantMassFit=massfit;} + + protected: + + Double_t fBeamEnergy; // Energy of the beam (required for the CS angle) + Double_t fInvMassFitLimitMin; // invariant mass spectrum fit limits + Double_t fInvMassFitLimitMax; + Double_t fPsiFitLimitMin; // psi fit limits + Double_t fPsiFitLimitMax; + Double_t fBckFitLimitMin; // bck fit limits + Double_t fBckFitLimitMax; + Bool_t fInvariantMassFit; // flag to perform or not inv. mass fit + + const char* fAnalysisType; //ESD or AOD based analysis + TList *fOutput; + + Float_t InvMass (Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, Float_t); + Float_t Rapidity (Float_t, Float_t); + Double_t CostCS (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t); + Double_t CostHE (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t); + void FitInvMass(TH1D* ); + + ClassDef(AliAnalysisTaskMuonDistributions,1); +}; + +#endif -- 2.39.3