From 9fe33e0bc7b8d42563d823ebf9f19c67be69dc1e Mon Sep 17 00:00:00 2001 From: martinez Date: Thu, 15 Apr 2010 21:58:05 +0000 Subject: [PATCH] New analysis taks for muon-dimuon analysis (Livio) --- PWG3/PWG3muonLinkDef.h | 1 + PWG3/libPWG3muon.pkg | 3 +- PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx | 635 +++++++++++++++++++ PWG3/muon/AliAnalysisTaskMuonTreeBuilder.h | 119 ++++ 4 files changed, 757 insertions(+), 1 deletion(-) create mode 100644 PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx create mode 100644 PWG3/muon/AliAnalysisTaskMuonTreeBuilder.h diff --git a/PWG3/PWG3muonLinkDef.h b/PWG3/PWG3muonLinkDef.h index 0394fd7cfc3..8c9933f7294 100644 --- a/PWG3/PWG3muonLinkDef.h +++ b/PWG3/PWG3muonLinkDef.h @@ -27,5 +27,6 @@ #pragma link C++ class AliMuonsHFHeader+; #pragma link C++ class AliAnalysisTaskSEMuonsHF+; #pragma link C++ class AliAnalysisTaskDimuonCFContainerBuilder+; +#pragma link C++ class AliAnalysisTaskMuonTreeBuilder+; #endif diff --git a/PWG3/libPWG3muon.pkg b/PWG3/libPWG3muon.pkg index 867e28da26d..595f4441a81 100644 --- a/PWG3/libPWG3muon.pkg +++ b/PWG3/libPWG3muon.pkg @@ -22,7 +22,8 @@ SRCS:= muon/AliAnalysisTaskESDMuonFilter.cxx \ muon/AliDimuInfoStoreMC.cxx \ muon/AliMuonsHFHeader.cxx \ muon/AliAnalysisTaskSEMuonsHF.cxx \ - muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx + muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx \ + muon/AliAnalysisTaskMuonTreeBuilder.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx b/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx new file mode 100644 index 00000000000..80928a04db1 --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx @@ -0,0 +1,635 @@ +#ifndef ALIANALYSISTASKMUONTREEBUILDER_CXX +#define ALIANALYSISTASKMUONTREEBUILDER_CXX + +#include "AliAnalysisTaskMuonTreeBuilder.h" +#include "AliMCEvent.h" +#include "AliESDMuonTrack.h" +#include "AliESDVertex.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliESDHeader.h" +#include "TParticle.h" +#include "TLorentzVector.h" +#include "TFile.h" +#include "TH1I.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "TChain.h" +#include "AliESDtrack.h" +#include "AliLog.h" +#include "AliESDtrack.h" +#include "AliESDInputHandler.h" +#include "TCanvas.h" +#include "AliPhysicsSelection.h" +#include "AliMultiplicity.h" + +// Analysis task for muon-dimuon analysis +// Works for real and MC events +// Works with the corresponding AddTask macro +// +// author: L. Bianchi - Universita' & INFN Torino + +ClassImp(AliAnalysisTaskMuonTreeBuilder) + +//__________________________________________________________________________ +AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder() : + fNevt(0), + fBeamEnergy(5000.), + fOutput(0x0), + fOutputTree(0x0), + fIsMC(kFALSE), + fIsSelected(kFALSE), + fNumMuonTracks(0), + fNumSPDTracklets(0), + fNumContributors(0), + fNumDimuons(0) +{ + // + //Default ctor + // +} +//___________________________________________________________________________ +AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const Char_t* name) : + AliAnalysisTaskSE(name), + fNevt(0), + fBeamEnergy(5000.), + fOutput(0x0), + fOutputTree(0x0), + fIsMC(kFALSE), + fIsSelected(kFALSE), + fNumMuonTracks(0), + fNumSPDTracklets(0), + fNumContributors(0), + fNumDimuons(0) +{ + // + // Constructor. Initialization of Inputs and Outputs + // + Info("AliAnalysisTaskMuonTreeBuilder","Calling Constructor"); + + DefineOutput(1,TTree::Class()); +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonTreeBuilder& AliAnalysisTaskMuonTreeBuilder::operator=(const AliAnalysisTaskMuonTreeBuilder& c) +{ + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c) ; + fNevt = c.fNevt ; + } + return *this; +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const AliAnalysisTaskMuonTreeBuilder& c) : + AliAnalysisTaskSE(c), + fNevt(c.fNevt), + fBeamEnergy(c.fBeamEnergy), + fOutput(c.fOutput), + fOutputTree(c.fOutputTree), + fIsMC(kFALSE), + fIsSelected(kFALSE), + fNumMuonTracks(c.fNumMuonTracks), + fNumSPDTracklets(c.fNumSPDTracklets), + fNumContributors(c.fNumContributors), + fNumDimuons(c.fNumDimuons) +{ + // + // Copy Constructor FIDUCIAL REGIONS? + // +} + +//___________________________________________________________________________ +AliAnalysisTaskMuonTreeBuilder::~AliAnalysisTaskMuonTreeBuilder() { + // + //destructor + // + Info("~AliAnalysisTaskMuonTreeBuilder","Calling Destructor"); +} + +//___________________________________________________________________________ +void AliAnalysisTaskMuonTreeBuilder::UserCreateOutputObjects(){ + +// +// Creating User-Defined Output Objects +// + + // TREE OUTPUT---------------------------------------------------------- + OpenFile(1); + fOutputTree = new TTree("krec","Tree of reconstructed muons"); + + fOutputTree->Branch("IsSelected",&fIsSelected,"IsSelected/O"); + fOutputTree->Branch("FiredTriggerClasses",fTrigClass,"FiredTriggerClasses/C"); + + fOutputTree->Branch("NumMuonTracks",&fNumMuonTracks,"NumMuonTracks/I"); + fOutputTree->Branch("NumContributors",&fNumContributors,"NumContributors/I"); + fOutputTree->Branch("NumSPDTracklets",&fNumSPDTracklets,"NumSPDTraclets/I"); + fOutputTree->Branch("Vertex",fVertex,"Vertex[3]/D"); + fOutputTree->Branch("pT",fpT,"pT[10]/D"); + fOutputTree->Branch("E",fE,"E[10]/D"); + fOutputTree->Branch("px",fpx,"px[10]/D"); + fOutputTree->Branch("py",fpy,"py[10]/D"); + fOutputTree->Branch("pz",fpz,"pz[10]/D"); + fOutputTree->Branch("pxUncorr",fpxUncorr,"pxUncorr[10]/D"); + fOutputTree->Branch("pyUncorr",fpyUncorr,"pyUncorr[10]/D"); + fOutputTree->Branch("pzUncorr",fpzUncorr,"pzUncorr[10]/D"); + fOutputTree->Branch("y",fy,"y[10]/D"); + fOutputTree->Branch("eta",feta,"eta[10]/D"); + fOutputTree->Branch("phi",fphi,"phi[10]/D"); + fOutputTree->Branch("MatchTrig",fMatchTrig,"MatchTrig[10]/I"); + fOutputTree->Branch("TrackChi2",fTrackChi2,"TrackChi2[10]/D"); + fOutputTree->Branch("MatchTrigChi2",fMatchTrigChi2,"MatchTrigChi2[10]/D"); + fOutputTree->Branch("DCA",fDCA,"DCA[10]/D"); + fOutputTree->Branch("Charge",fCharge,"Charge[10]/S"); + fOutputTree->Branch("MuFamily",fMuFamily,"MuFamily[10]/I"); + fOutputTree->Branch("RAtAbsEnd",fRAtAbsEnd,"RAtAbsEnd[10]/D"); + + fOutputTree->Branch("NumDimuons",&fNumDimuons,"NumDimuons/I"); + fOutputTree->Branch("DimuConstituent",fDimuonConstituent,"DimuonConstituent[45][2]/I"); + fOutputTree->Branch("pTdimuon",fpTdimuon,"pTdimuon[45]/D"); + fOutputTree->Branch("pxdimuon",fpxdimuon,"pxdimuon[45]/D"); + fOutputTree->Branch("pydimuon",fpydimuon,"pydimuon[45]/D"); + fOutputTree->Branch("pzdimuon",fpzdimuon,"pzdimuon[45]/D"); + fOutputTree->Branch("ydimuon",fydimuon,"ydimuon[45]/D"); + fOutputTree->Branch("Imassdimuon",fiMassdimuon,"iMassdimuon[45]/D"); + fOutputTree->Branch("costCS",fcostCS,"costCS[45]/D"); + fOutputTree->Branch("costHE",fcostHE,"costHE[45]/D"); + fOutputTree->Branch("phiCS",fphiCS,"phiCS[45]/D"); + fOutputTree->Branch("phiHE",fphiHE,"phiHE[45]/D"); + + fOutputTree->Branch("PDG",fPDG,"PDG[10]/I"); + fOutputTree->Branch("PDGmother",fPDGmother,"PDGmother[10]/I"); + fOutputTree->Branch("PDGdimu",fPDGdimu,"PDGdimu[45]/I"); + +} + + + +//_________________________________________________ +void AliAnalysisTaskMuonTreeBuilder::UserExec(Option_t *) +{ +// +// User Exec +// + + fNevt++; + + + fNumMuonTracks=0; + fNumSPDTracklets=666.; + fNumContributors=666.; + fNumDimuons=0; + fIsSelected=kFALSE; + fVertex[0]=666.; fVertex[1]=666.; fVertex[2]=666.; + for(Int_t i=0; i<10;i++){ + fpT[i]=666.; + fE[i]=666.; + fpx[i]=666; + fpy[i]=666; + fpz[i]=666; + fpxUncorr[i]=666; + fpyUncorr[i]=666; + fpzUncorr[i]=666; + fy[i]=666.; + feta[i]=666.; + fphi[i]=666.; + fMatchTrig[i]=666.; + fTrackChi2[i]=666.; + fMatchTrigChi2[i]=666.; + fDCA[i]=666.; + fPDG[i]=666; + fPDGmother[i]=666; + fCharge[i]=666; + fMuFamily[i]=666; + fRAtAbsEnd[i]=666; + } + for(Int_t i=0; i<45;i++){ + fDimuonConstituent[i][0]=666; fDimuonConstituent[i][1]=666; + fpTdimuon[i]=666.; + fpxdimuon[i]=666.; + fpydimuon[i]=666.; + fpzdimuon[i]=666.; + fydimuon[i]=666.; + fiMassdimuon[i]=666.; + fcostCS[i]=666.; + fcostHE[i]=666.; + fphiCS[i]=666.; + fphiHE[i]=666.; + fPDGdimu[i]=666; + } + + +//////// +//// ESD +//////// + + AliESDEvent *fESD = 0x0; + AliMCEvent* mcEvent = 0x0; + + if(fIsMC){ + if (!fMCEvent) { + Error("UserExec","NO MC EVENT FOUND!"); + return; + } + } + + fESD = dynamic_cast(InputEvent()); + + fIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); + + fNumMuonTracks = fESD->GetNumberOfMuonTracks() ; + Int_t loopEnd = fNumMuonTracks; + if(!fIsMC) { + TString cla = fESD->GetFiredTriggerClasses(); + sprintf(fTrigClass,"%s",cla.Data()); + } + + if(fNumMuonTracks>0 && fIsMC){ + mcEvent = MCEvent(); + } + fNumContributors = fESD->GetPrimaryVertexSPD()->GetNContributors(); + fNumSPDTracklets = fESD->GetMultiplicity()->GetNumberOfTracklets(); + fVertex[0]=fESD->GetPrimaryVertexSPD()->GetX(); + fVertex[1]=fESD->GetPrimaryVertexSPD()->GetY(); + fVertex[2]=fESD->GetPrimaryVertexSPD()->GetZ(); + printf("fVertex : %f - %f - %f\n",fVertex[0],fVertex[1],fVertex[2]); + + Int_t numdimu = 0; + for (Int_t j = 0; jGetMuonTrack(j))); + if (!mu1->ContainTrackerData()) {fNumMuonTracks=fNumMuonTracks-1; continue;} + Double_t charge1 = mu1->Charge(); + fCharge[j] = mu1->Charge(); + fpT[j] = mu1->Pt(); + fpx[j] = mu1->Px(); + fpy[j] = mu1->Py(); + fpz[j] = mu1->Pz(); + fpxUncorr[j] = mu1->PxUncorrected(); + fpyUncorr[j] = mu1->PyUncorrected(); + fpzUncorr[j] = mu1->PzUncorrected(); + fy[j] = mu1->Y(); + feta[j] = mu1->Eta(); + fphi[j] = Phideg(mu1->Phi()); + Double_t emu1 = mu1->E(); + fE[j] = emu1; + fDCA[j] = mu1->GetDCA(); + fTrackChi2[j] = mu1->GetChi2(); + fMatchTrig[j] = mu1->GetMatchTrigger(); + fMatchTrigChi2[j]= mu1->GetChi2MatchTrigger(); + fRAtAbsEnd[j]=mu1->GetRAtAbsorberEnd(); + + AliMCParticle* mcTrack = 0x0; + if(fIsMC){ + if(mu1->GetLabel()==-1) continue; + mcTrack = (AliMCParticle*)mcEvent->GetTrack(mu1->GetLabel()); + fPDG[j] = mcTrack->PdgCode(); + if (mcTrack->GetMother()==-1) continue; + fPDGmother[j] = ((AliMCParticle*)mcEvent->GetTrack(mcTrack->GetMother()))->PdgCode(); + if (TMath::Abs(fPDG[j])==13) fMuFamily[j] = FindMuFamily(mcTrack,mcEvent); + } + for (Int_t jj = j+1; jjGetMuonTrack(jj))); + if (!mu2->ContainTrackerData()) continue; + + Double_t pxmu2 = mu2->Px(); + Double_t pymu2 = mu2->Py(); + Double_t pzmu2 = mu2->Pz(); + Double_t emu2 = mu2->E(); + //Double_t charge2= mu2->Charge(); + + fpTdimuon[numdimu] = TMath::Sqrt((fpx[j]+pxmu2)*(fpx[j]+pxmu2)+(fpy[j]+pymu2)*(fpy[j]+pymu2)); + fpxdimuon[numdimu] = fpx[j]+pxmu2; + fpydimuon[numdimu] = fpy[j]+pymu2; + fpzdimuon[numdimu] = fpz[j]+pzmu2; + fydimuon[numdimu] = Rap((emu1+emu2),(fpz[j]+pzmu2)); + fiMassdimuon[numdimu] = Imass(emu1,fpx[j],fpy[j],fpz[j],emu2,pxmu2,pymu2,pzmu2); + fcostCS[numdimu]=CostCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2); + fcostHE[numdimu]=CostHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2); + fphiCS[numdimu] = PhiCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2); + fphiHE[numdimu] = PhiHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2); + + fDimuonConstituent[numdimu][0]=j; fDimuonConstituent[numdimu][1]=jj; + + numdimu++; + + if(fIsMC){ + if(mu2->GetLabel()==-1) continue; + if(TMath::Abs(fPDG[j])==13 && TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()))->PdgCode())==13) fPDGdimu[numdimu-1]=FindDimuFamily(mcTrack,(AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()),mcEvent); + else fPDGdimu[numdimu-1]=-3; + } + + delete mu2; + } + fNumDimuons=numdimu; + + + delete mu1; + } + fOutputTree->Fill(); + PostData(1,fOutputTree); +} + +//________________________________________________________________________ +Int_t AliAnalysisTaskMuonTreeBuilder::FindDimuFamily(AliMCParticle* mcTrack1,AliMCParticle* mcTrack2, AliMCEvent* mcEvent) const +{ + // finds the family of the dimuon (works only if the 2 muons are real muons (not hadrons)) + Int_t familynumber; + + if(mcTrack1->GetMother()==mcTrack2->GetMother()) familynumber = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mcTrack1->GetMother()))->PdgCode()); + else{ + Int_t familymu1 = FindMuFamily(mcTrack1,mcEvent); + Int_t familymu2 = FindMuFamily(mcTrack2,mcEvent); + if(familymu1==5 && familymu2==5) familynumber=5; //bb dimuon + else if(familymu1==4 && familymu2==4) familynumber=4; //cc dimuon + else if((familymu1==4 && familymu2==5)||(familymu2==4 && familymu1==5)) familynumber=45; //bc dimuon + else if (familymu1==-2 || familymu2==-2 || familymu1==-3 || familymu2==-3) familynumber=-2; //hadron dimuon (at least 1 hadron involved) + else familynumber=-1; + } + + return familynumber; +} + +//________________________________________________________________________ +Int_t AliAnalysisTaskMuonTreeBuilder::FindMuFamily(AliMCParticle* mcTrack, AliMCEvent* mcEvent) const +{ + // finds the family of the muon + Int_t imother = mcTrack->GetMother(); + if ( imother<0 ) return -1; // Drell-Yan Muon + + Int_t igrandma = imother; + + AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother); + Int_t motherPdg = motherPart->PdgCode(); + + // Track is an heavy flavor muon + Int_t absPdg = TMath::Abs(motherPdg); + if(absPdg/100==5 || absPdg/1000==5) return 5; + if(absPdg/100==4 || absPdg/1000==4){ + Int_t newMother = -1; + igrandma = imother; + Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode()); + while ( absGrandMotherPdg > 10 ) { + igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother(); + if( igrandma < 0 ) break; + absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode()); + } + + if (absGrandMotherPdg==5) newMother = 5; // Charm from beauty + else if (absGrandMotherPdg==4) newMother = 4; + + if(newMother<0) { + //AliWarning("Mother not correctly found! Set to charm!\n"); + newMother = 4; + } + + return newMother; + } + + Int_t nPrimaries = mcEvent->Stack()->GetNprimary(); + // Track is a bkg. muon + if (imother0) {cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());} + else {cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());} + return cost; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskMuonTreeBuilder::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) +{ +// calculate costHE + 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; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskMuonTreeBuilder::PhiCS(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) +{ +// calculate phiCS + TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,yaxisCS, xaxisCS, 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(); + yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit(); + xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit(); + + Double_t phi=0.; + if(charge1>0) { + phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS))); + } else { + phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS))); + } + if (phi>TMath::Pi()) phi=phi-TMath::Pi(); + return phi; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskMuonTreeBuilder::PhiHE(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) +{ +// calculate phiHE + TLorentzVector PMu1Lab, PMu2Lab, PProjLab, PTargLab, PDimuLab; // In the lab. frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,xaxis, yaxis,zaxis; + Double_t mp=0.93827231; + + // + // --- Get the muons parameters in the LAB frame + // + PMu1Lab.SetPxPyPzE(px1,py1,pz1,e1); + PMu2Lab.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the LAB frame + // + PDimuLab=PMu1Lab+PMu2Lab; + + zaxis=(PDimuLab.Vect()).Unit(); + // + // --- Translate the muon parameters in the dimuon rest frame + // + beta=(-1./PDimuLab.E())*PDimuLab.Vect(); + + PProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile + PTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio + + PProjDimu=PProjLab; + PTargDimu=PTargLab; + + PProjDimu.Boost(beta); + PTargDimu.Boost(beta); + + yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit(); + xaxis=(yaxis.Cross(zaxis)).Unit(); + + PMu1Dimu=PMu1Lab; + PMu2Dimu=PMu2Lab; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + + Double_t phi=0.; + if(charge1 > 0) { + phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis)); + } else { + phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis)); + } + return phi; +} + +//________________________________________________________________________ +void AliAnalysisTaskMuonTreeBuilder::Terminate(Option_t *) +{ +// Terminate +} + +#endif diff --git a/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.h b/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.h new file mode 100644 index 00000000000..a1548e0f9a7 --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonTreeBuilder.h @@ -0,0 +1,119 @@ +#ifndef ALIANALYSISTASKMUONTREEBUILDER_H +#define ALIANALYSISTASKMUONTREEBUILDER_H + +#include "AliAnalysisTaskSE.h" +#include "TMath.h" + +// Analysis task for muon-dimuon analysis +// Works for real and MC events +// author: L. Bianchi - Universita' & INFN Torino + +class TH1I; +class TParticle ; +class TLorentzVector ; +//class TVector3 ; +class TFile ; +class AliStack ; +class AliESDtrack; +class AliVParticle; +class AliMCParticle; +class AliMCEvent; + + +class AliAnalysisTaskMuonTreeBuilder : public AliAnalysisTaskSE { + public: + + AliAnalysisTaskMuonTreeBuilder(); + AliAnalysisTaskMuonTreeBuilder(const Char_t* name); + AliAnalysisTaskMuonTreeBuilder& operator= (const AliAnalysisTaskMuonTreeBuilder& c); + AliAnalysisTaskMuonTreeBuilder(const AliAnalysisTaskMuonTreeBuilder& c); + virtual ~AliAnalysisTaskMuonTreeBuilder(); + + // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects + void UserExec(Option_t *option); + void Terminate(Option_t *); + void UserCreateOutputObjects(); + + + // Data types + void SetIsMC (Bool_t flagMC) {fIsMC=flagMC;} + void SetBeamEnergy (Double_t en) {fBeamEnergy=en;} + + protected: + + Double_t fNevt ; // event counter + Double_t fBeamEnergy ; // Energy of the beam (required for the CS angle) + TList *fOutput ; // output + TTree *fOutputTree ; //! tree output + + Bool_t fIsMC; + + Bool_t fIsSelected; + char fTrigClass[100]; + + Int_t fNumMuonTracks ; // variables for single mu + Int_t fNumSPDTracklets; + Int_t fNumContributors; + Double_t fVertex[3]; + Double_t fpT[10]; + Double_t fE[10]; + Double_t fpx[10]; + Double_t fpy[10]; + Double_t fpz[10]; + Double_t fpxUncorr[10]; + Double_t fpyUncorr[10]; + Double_t fpzUncorr[10]; + Double_t fy[10]; + Double_t feta[10]; + Double_t fphi[10]; + Int_t fMatchTrig[10]; + Double_t fTrackChi2[10]; + Double_t fMatchTrigChi2[10]; + Double_t fDCA[10]; + Short_t fCharge[10]; + Int_t fMuFamily[10]; + Double_t fRAtAbsEnd[10]; + + Int_t fNumDimuons; // variables for dimuons + Int_t fDimuonConstituent[45][2]; + Double_t fpTdimuon[45]; + Double_t fpxdimuon[45]; + Double_t fpydimuon[45]; + Double_t fpzdimuon[45]; + Double_t fydimuon[45]; + Double_t fiMassdimuon[45]; + Double_t fcostCS[45]; + Double_t fcostHE[45]; + Double_t fphiCS[45]; + Double_t fphiHE[45]; + + //Int_t fIsPrimary[10]; + Int_t fPDG[10]; + Int_t fPDGmother[10]; + Int_t fPDGdimu[45]; + + + + + + Int_t FindDimuFamily(AliMCParticle* mcTrack1,AliMCParticle* mcTrack2, AliMCEvent* mcEvent) const; + Int_t FindMuFamily(AliMCParticle* mcTrack, AliMCEvent* mcEvent) const; + Double_t Imass (Double_t e1, Double_t px1, Double_t py1, Double_t pz1, Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const; + Double_t Rap (Double_t e, Double_t pz) const; +// Double_t Imass(Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t) const; +// Double_t Rap(Double_t, Double_t) const; + Double_t Phideg(Double_t phi) const; + +// 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); +// Double_t PhiCS (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t); +// Double_t PhiHE (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t); + Double_t 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); + Double_t 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); + Double_t PhiCS (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); + Double_t PhiHE (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); + + ClassDef(AliAnalysisTaskMuonTreeBuilder,1); +}; + +#endif -- 2.43.0