X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FMUONmassPlot_ESD.C;h=3acd495f518cfbe791ba042daffa032bebd0b746;hb=4eef97dc01d3ee6f5ff6971751592dc96b8b5ccb;hp=21af0a13e271431b1924787b30a7e6eddd672ac3;hpb=f57d136a5868c5d3576e89d0594e729ae59566c2;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/MUONmassPlot_ESD.C b/MUON/MUONmassPlot_ESD.C index 21af0a13e27..3acd495f518 100644 --- a/MUON/MUONmassPlot_ESD.C +++ b/MUON/MUONmassPlot_ESD.C @@ -10,54 +10,59 @@ #include "TParticle.h" #include "TTree.h" #include +#include +#include // STEER includes -#include "AliRun.h" -#include "AliRunLoader.h" -#include "AliHeader.h" -#include "AliLoader.h" -#include "AliStack.h" +#include "AliLog.h" #include "AliMagF.h" -#include "AliESD.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliTracker.h" +#include "AliESDMuonTrack.h" // MUON includes -#include "AliESDMuonTrack.h" +#include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" +#include "AliMUONESDInterface.h" #endif -// -// Macro MUONmassPlot.C for ESD -// Ch. Finck, Subatech, April. 2004 -// - -// macro to make invariant mass plots -// for combinations of 2 muons with opposite charges, -// from root file "MUON.tracks.root" containing the result of track reconstruction. -// Histograms are stored on the "MUONmassPlot.root" file. -// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...) -// using Invariant Mass for rapidity. - -// Arguments: -// FirstEvent (default 0) -// LastEvent (default 0) -// ResType (default 553) -// 553 for Upsilon, anything else for J/Psi -// Chi2Cut (default 100) -// to keep only tracks with chi2 per d.o.f. < Chi2Cut -// PtCutMin (default 1) -// to keep only tracks with transverse momentum > PtCutMin -// PtCutMax (default 10000) -// to keep only tracks with transverse momentum < PtCutMax -// massMin (default 9.17 for Upsilon) -// & massMax (default 9.77 for Upsilon) -// to calculate the reconstruction efficiency for resonances with invariant mass -// massMin < mass < massMax. - -// Add parameters and histograms for analysis - -Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 10000, - char* esdFileName = "AliESDs.root", Int_t ResType = 553, + +/// \ingroup macros +/// \file MUONmassPlot_ESD.C +/// \brief Macro MUONefficiency.C for ESD +/// +/// \author Ch. Finck, Subatech, April. 2004 +/// +/// +/// Macro to make invariant mass plots +/// for combinations of 2 muons with opposite charges, +/// from root file "MUON.tracks.root" containing the result of track reconstruction. +/// Histograms are stored on the "MUONmassPlot.root" file. +/// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...) +/// using Invariant Mass for rapidity. +/// +/// Add parameters and histograms for analysis + +Bool_t MUONmassPlot(Int_t ExtrapToVertex = -1, char* geoFilename = "geometry.root", + Int_t FirstEvent = 0, Int_t LastEvent = 10000, char* esdFileName = "AliESDs.root", Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000., Float_t massMin = 9.17,Float_t massMax = 9.77) { +/// \param ExtrapToVertex (default -1) +/// - <0: no extrapolation; +/// - =0: extrapolation to (0,0,0); +/// - >0: extrapolation to ESDVertex if available, else to (0,0,0) +/// \param FirstEvent (default 0) +/// \param LastEvent (default 0) +/// \param ResType 553 for Upsilon, anything else for J/Psi (default 553) +/// \param Chi2Cut to keep only tracks with chi2 per d.o.f. < Chi2Cut (default 100) +/// \param PtCutMin to keep only tracks with transverse momentum > PtCutMin (default 1) +/// \param PtCutMax to keep only tracks with transverse momentum < PtCutMax (default 10000) +/// \param massMin (default 9.17 for Upsilon) +/// \param massMax (default 9.77 for Upsilon); +/// to calculate the reconstruction efficiency for resonances with invariant mass +/// massMin < mass < massMax. + cout << "MUONmassPlot " << endl; cout << "FirstEvent " << FirstEvent << endl; cout << "LastEvent " << LastEvent << endl; @@ -84,7 +89,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.); TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.); TH1F *hInvMassRes; - TH1F *hPrimaryVertex = new TH1F("hPrimaryVertex","SPD reconstructed Z vertex",120,-12,12); + TH1F *hPrimaryVertex = new TH1F("hPrimaryVertex","SPD reconstructed Z vertex",150,-15,15); if (ResType == 553) { hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.); @@ -109,34 +114,39 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // Float_t UpsilonMass = 9.46037; // Float_t JPsiMass = 3.097; - Double_t thetaX, thetaY, pYZ; + Int_t fCharge1, fCharge2; Double_t fPxRec1, fPyRec1, fPzRec1, fE1; Double_t fPxRec2, fPyRec2, fPzRec2, fE2; - Int_t fCharge, fCharge2; Int_t ntrackhits, nevents; Double_t fitfmin; - Double_t fZVertex; - + Double_t fZVertex=0; + Double_t fYVertex=0; + Double_t fXVertex=0; + Double_t errXVtx=0; + Double_t errYVtx=0; TLorentzVector fV1, fV2, fVtot; - // set off mag field - AliMagF::SetReadField(kFALSE); - - // open run loader and load gAlice, kinematics and header - AliRunLoader* runLoader = AliRunLoader::Open(filename); - if (!runLoader) { - Error("MUONmass_ESD", "getting run loader from file %s failed", - filename); - return kFALSE; - } - - if (!gAlice) { - Error("MUONmass_ESD", "no galice object found"); - return kFALSE; + // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex) + if (!gGeoManager) { + TGeoManager::Import(geoFilename); + if (!gGeoManager) { + Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename); + return kFALSE; + } } + // set mag field + // waiting for mag field in CDB + if (!TGeoGlobalMagField::Instance()->GetField()) { + printf("Loading field map...\n"); + AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG); + TGeoGlobalMagField::Instance()->SetField(field); + } + // set the magnetic field for track extrapolations + AliMUONTrackExtrap::SetField(); + // open the ESD file TFile* esdFile = TFile::Open(esdFileName); @@ -145,26 +155,22 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t return kFALSE; } - AliESD* esd = new AliESD(); + AliESDEvent* esd = new AliESDEvent(); TTree* tree = (TTree*) esdFile->Get("esdTree"); if (!tree) { Error("CheckESD", "no ESD tree found"); return kFALSE; } - tree->SetBranchAddress("ESD", &esd); +// tree->SetBranchAddress("ESD", &esd); + esd->ReadFromTree(tree); + nevents = (Int_t)tree->GetEntries(); - AliESDVertex* Vertex = (AliESDVertex*) esd->AliESD::GetVertex(); + AliMUONTrackParam trackParam; - runLoader->LoadHeader(); - nevents = runLoader->GetNumberOfEvents(); - // Loop over events for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) { - // get current event - runLoader->GetEvent(iEvent); - // get the event summary data tree->GetEvent(iEvent); if (!esd) { @@ -173,7 +179,14 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t } // get the SPD reconstructed vertex (vertexer) and fill the histogram - fZVertex = Vertex->GetZv(); + AliESDVertex* Vertex = (AliESDVertex*) esd->GetVertex(); + if (Vertex->GetNContributors()) { + fZVertex = Vertex->GetZv(); + fYVertex = Vertex->GetYv(); + fXVertex = Vertex->GetXv(); + errXVtx = Vertex->GetXRes(); + errYVtx = Vertex->GetYRes(); + } hPrimaryVertex->Fill(fZVertex); Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; @@ -184,20 +197,26 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // loop over all reconstructed tracks (also first track of combination) for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); - - thetaX = muonTrack->GetThetaX(); - thetaY = muonTrack->GetThetaY(); - - pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum()); - fPzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY)); - fPxRec1 = fPzRec1 * TMath::Tan(thetaX); - fPyRec1 = fPzRec1 * TMath::Tan(thetaY); - fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); - - fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1); - fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1); - + // skip ghosts + if (!esd->GetMuonTrack(iTrack)->ContainTrackerData()) continue; + + AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack))); + + // extrapolate to vertex if required and available + if (ExtrapToVertex > 0 && Vertex->GetNContributors()) { + AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack, trackParam); + AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx); + AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack); // put the new parameters in this copy of AliESDMuonTrack + } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){ + AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack, trackParam); + AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.); + AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack); // put the new parameters in this copy of AliESDMuonTrack + } + + fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); + + muonTrack->LorentzP(fV1); + ntrackhits = muonTrack->GetNHit(); fitfmin = muonTrack->GetChi2(); @@ -213,7 +232,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // chi2 per d.o.f. Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5); // printf(" px %f py %f pz %f NHits %d Norm.chi2 %f charge %d\n", -// fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge); +// fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge1); // condition for good track (Chi2Cut and PtCut) @@ -224,32 +243,38 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t hPMuon->Fill(p1); hChi2PerDof->Fill(ch1); hRapMuon->Fill(rapMuon1); - if (fCharge > 0) { + if (fCharge1 > 0) { hPtMuonPlus->Fill(pt1); - hThetaPhiPlus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415); + hThetaPhiPlus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi()); } else { hPtMuonMinus->Fill(pt1); - hThetaPhiMinus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415); + hThetaPhiMinus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi()); } // loop over second track of combination for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) { - AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2); - - thetaX = muonTrack->GetThetaX(); - thetaY = muonTrack->GetThetaY(); - - pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum()); - fPzRec2 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY)); - fPxRec2 = fPzRec2 * TMath::Tan(thetaX); - fPyRec2 = fPzRec2 * TMath::Tan(thetaY); - fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); + // skip ghosts + if (!esd->GetMuonTrack(iTrack2)->ContainTrackerData()) continue; + + AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2))); + + // extrapolate to vertex if required and available + if (ExtrapToVertex > 0 && Vertex->GetNContributors()) { + AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack2, trackParam); + AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx); + AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack2); // put the new parameters in this copy of AliESDMuonTrack + } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){ + AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack2, trackParam); + AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.); + AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack2); // put the new parameters in this copy of AliESDMuonTrack + } + + fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum())); - fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2); - fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2); + muonTrack2->LorentzP(fV2); - ntrackhits = muonTrack->GetNHit(); - fitfmin = muonTrack->GetChi2(); + ntrackhits = muonTrack2->GetNHit(); + fitfmin = muonTrack2->GetChi2(); // transverse momentum Float_t pt2 = fV2.Pt(); @@ -261,7 +286,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) { // condition for opposite charges - if ((fCharge * fCharge2) == -1) { + if ((fCharge1 * fCharge2) == -1) { // invariant mass fVtot = fV1 + fV2; @@ -273,24 +298,26 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt()); Int_t ptTrig; if (ResType == 553) - ptTrig = 0x400;// mask for Hpt unlike sign pair + ptTrig = 0x20;// mask for Hpt unlike sign pair else - ptTrig = 0x200;// mask for Lpt unlike sign pair + ptTrig = 0x10;// mask for Lpt unlike sign pair if (esd->GetTriggerMask() & ptTrig) NbTrigger++; if (invMass > massMin && invMass < massMax) { EventInMass++; - if (muonTrack->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger + if (muonTrack2->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger EventInMassMatch++; hRapResonance->Fill(fVtot.Rapidity()); hPtResonance->Fill(fVtot.Pt()); } - } // if (fCharge * fCharge2) == -1) + } // if (fCharge1 * fCharge2) == -1) } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) + delete muonTrack2; } // for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++) } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) ) + delete muonTrack; } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++) hNumberOfTrack->Fill(nTracks);