X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FMUONefficiency.C;h=e799880cd7cb14e8a88d603ecf95f856085d69fa;hb=96c0c85603eba05064784139723787825bf1c300;hp=204f83efd5e0ae8191da5ddbba21db917bea10fc;hpb=3f6289e9ae207cbf82d0e94e56e6050fcab90ae3;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/MUONefficiency.C b/MUON/MUONefficiency.C index 204f83efd5e..e799880cd7c 100644 --- a/MUON/MUONefficiency.C +++ b/MUON/MUONefficiency.C @@ -15,61 +15,72 @@ /* $Id$ */ -// Macro (upgraded version of MUONmassPlot_ESD.C, better handling of Jpsi) to make : -// 1) Ntuple (Ktuple) containing Upsilon kinematics variables (from kinematics.root files) -// 2) Ntuple (ESDtuple) containing Upsilon kinematics variables from reconstruction and -// combinations of 2 muons with opposite charges (ESDtupleBck will be used later) -// 3) Some QA histograms -// Ntuple are stored in the file MUONefficiency.root and ESD tree and QA histograms in AliESDs.root +/// \ingroup macros +/// \file MUONefficiency.C +/// \brief add brief description +/// +/// Macro (upgraded version of MUONmassPlot_ESD.C, better handling of Jpsi) to make : +/// - Ntuple (Ktuple) containing Upsilon kinematics variables (from kinematics.root files) +/// - Ntuple (ESDtuple) containing Upsilon kinematics variables from reconstruction and +/// combinations of 2 muons with opposite charges (ESDtupleBck will be used later) +/// - Some QA histograms +/// Ntuple are stored in the file MUONefficiency.root and ESD tree and QA histograms in AliESDs.root +/// +/// \author Christophe Suire, IPN Orsay -// Christophe Suire, IPN Orsay +#if !defined(__CINT__) || defined(__MAKECINT__) -// Arguments: -// FirstEvent (default 0) -// LastEvent (default 1.e6) -// ResType (default 553) -// 553 for Upsilon, 443 for J/Psi -// Chi2Cut (default 100) -// to keep only tracks with chi2 per d.o.f. < Chi2Cut - +// MUON includes +#include "AliMUONCDB.h" +#include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" +#include "AliMUONESDInterface.h" +// STEER includes +#include "AliRun.h" +#include "AliRunLoader.h" +#include "AliHeader.h" +#include "AliLoader.h" +#include "AliStack.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliCDBManager.h" +#include "AliESDMuonTrack.h" -#if !defined(__CINT__) || defined(__MAKECINT__) // ROOT includes #include "TTree.h" #include "TNtuple.h" -#include "TBranch.h" -#include "TClonesArray.h" #include "TLorentzVector.h" #include "TFile.h" #include "TH1.h" #include "TH2.h" #include "TParticle.h" -#include "TTree.h" #include "TString.h" #include +#include +#include +#include +#include -// STEER includes -#include "AliRun.h" -#include "AliRunLoader.h" -#include "AliHeader.h" -#include "AliLoader.h" -#include "AliStack.h" -#include "AliMagFMaps.h" -#include "AliESD.h" -#include "AliTracker.h" - -// MUON includes -#include "AliMUONTrackParam.h" -#include "AliESDMuonTrack.h" #endif +Bool_t MUONefficiency(char* filename = "generated/galice.root", char* esdFileName = "AliESDs.root", + char* geoFilename = "geometry.root", char* ocdbPath = "local://$ALICE_ROOT/OCDB", + Int_t ExtrapToVertex = -1, Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000 ) +{ +/// \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 ResType 553 for Upsilon, 443 for J/Psi (default 553) +/// \param FirstEvent (default 0) +/// \param LastEvent (default 1.e6) +/// \param Chi2Cut to keep only tracks with chi2 per d.o.f. < Chi2Cut (default 100) -Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000, - char* esdFileName = "AliESDs.root", char* filename = "galice.root") -{ // MUONefficiency starts + + // MUONefficiency starts Double_t MUON_MASS = 0.105658369; Double_t UPSILON_MASS = 9.4603 ; @@ -164,8 +175,9 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven Double_t fXVertex=0; Double_t fYVertex=0; Double_t fZVertex=0; + Double_t errXVtx=0; + Double_t errYVtx=0; - Double_t thetaX, thetaY, pYZ; Double_t fPxRec1, fPyRec1, fPzRec1, fE1; Double_t fPxRec2, fPyRec2, fPzRec2, fE2; Int_t fCharge1, fCharge2; @@ -176,11 +188,44 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven TLorentzVector fV1, fV2, fVtot; - // set mag field - // waiting for mag field in CDB - printf("Loading field map...\n"); - AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k4kG); - AliTracker::SetFieldMap(field, kFALSE); + // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex) + if (!gGeoManager) { + TGeoManager::Import(geoFilename); + if (!gGeoManager) { + Error("MUONefficiency", "getting geometry from file %s failed", geoFilename); + return kFALSE; + } + } + + // open the ESD file + TFile* esdFile = TFile::Open(esdFileName); + if (!esdFile || !esdFile->IsOpen()) { + Error("MUONefficiency", "opening ESD file %s failed", esdFileName); + return kFALSE; + } + + AliESDEvent* esd = new AliESDEvent(); + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if (!tree) { + Error("MUONefficiency", "no ESD tree found"); + return kFALSE; + } + esd->ReadFromTree(tree); + + // get run number + if (tree->GetEvent(0) <= 0) { + Error("MUONefficiency", "no ESD object found for event 0"); + return kFALSE; + } + Int_t runNumber = esd->GetRunNumber(); + + // load necessary data from OCDB + AliCDBManager::Instance()->SetDefaultStorage(ocdbPath); + AliCDBManager::Instance()->SetRun(runNumber); + if (!AliMUONCDB::LoadField()) return kFALSE; + + // set the magnetic field for track extrapolations + AliMUONTrackExtrap::SetField(); // open run loader and load gAlice, kinematics and header AliRunLoader* runLoader = AliRunLoader::Open(filename); @@ -195,23 +240,13 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven Error("MUONefficiency", "no galice object found"); return kFALSE; } - - // open the ESD file - TFile* esdFile = TFile::Open(esdFileName); - if (!esdFile || !esdFile->IsOpen()) { - Error("MUONefficiency", "opening ESD file %s failed", esdFileName); - return kFALSE; - } - - AliESD* esd = new AliESD(); - TTree* tree = (TTree*) esdFile->Get("esdTree"); - if (!tree) { - Error("CheckESD", "no ESD tree found"); + + runLoader->LoadHeader(); + if (runNumber != runLoader->GetHeader()->GetRun()) { + Error("MUONefficiency", "mismatch between run number from ESD and from runLoader"); return kFALSE; } - tree->SetBranchAddress("ESD", &esd); - runLoader->LoadHeader(); nevents = runLoader->GetNumberOfEvents(); AliMUONTrackParam trackParam; @@ -291,18 +326,19 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven // get the event summary data - tree->GetEvent(iEvent); - if (!esd) { + if (tree->GetEvent(iEvent) <= 0) { Error("CheckESD", "no ESD object found for event %d", iEvent); return kFALSE; } // get the SPD reconstructed vertex (vertexer) and fill the histogram - AliESDVertex* Vertex = (AliESDVertex*) esd->AliESD::GetVertex(); - if (Vertex) { + 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); @@ -317,13 +353,20 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven // loop over all reconstructed tracks (also first track of combination) for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); + // skip ghosts + if (!esd->GetMuonTrack(iTrack)->ContainTrackerData()) continue; + + AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack))); - if (!Vertex) { - //re-extrapolate to vertex, if not kown before. - trackParam.GetParamFrom(*muonTrack); - trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex); - trackParam.SetParamFor(*muonTrack); + // 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 } // Trigger @@ -334,17 +377,9 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven else track1TriggerChi2 = 0. ; - 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); fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); - fE1 = TMath::Sqrt(MUON_MASS * MUON_MASS + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1); - fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1); + muonTrack->LorentzP(fV1); ntrackhits = muonTrack->GetNHit(); fitfmin = muonTrack->GetChi2(); @@ -361,7 +396,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven // chi2 per d.o.f. Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5); - if (PRINTLEVEL > 5 ) printf(" px %f py %f pz %f pt %f NHits %d Norm.chi2 %f charge %d\n",fPxRec1, fPyRec1, fPzRec1, pt1, ntrackhits, ch1, fCharge1); + if (PRINTLEVEL > 5 ) printf(" px %f py %f pz %f pt %f NHits %d Norm.chi2 %f charge %d\n",fV1.Px(), fV1.Py(), fV1.Pz(), pt1, ntrackhits, ch1, fCharge1); if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) { // condition for good track (Chi2Cut and PtCut) @@ -375,43 +410,43 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven 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); + // skip ghosts + if (!esd->GetMuonTrack(iTrack2)->ContainTrackerData()) continue; + + AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2))); - if (!Vertex) { - trackParam.GetParamFrom(*muonTrack); - trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex); - trackParam.SetParamFor(*muonTrack); + // 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 } - track2Trigger = muonTrack->GetMatchTrigger(); + track2Trigger = muonTrack2->GetMatchTrigger(); if (track2Trigger) - track2TriggerChi2 = muonTrack->GetChi2MatchTrigger(); + track2TriggerChi2 = muonTrack2->GetChi2MatchTrigger(); else track2TriggerChi2 = 0. ; - thetaX = muonTrack->GetThetaX(); - thetaY = muonTrack->GetThetaY(); + fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum())); - 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())); + muonTrack2->LorentzP(fV2); - fE2 = TMath::Sqrt(MUON_MASS * MUON_MASS + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2); - fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2); - - ntrackhits = muonTrack->GetNHit(); - fitfmin = muonTrack->GetChi2(); + ntrackhits = muonTrack2->GetNHit(); + fitfmin = muonTrack2->GetChi2(); // transverse momentum Float_t pt2 = fV2.Pt(); @@ -448,11 +483,9 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven //trigger info if (ResType == 553) - ptTrig = 0x400;// mask for Hpt unlike sign pair + ptTrig = 0x08;// mask for Hpt unlike sign pair else if (ResType == 443) - ptTrig = 0x800;// mask for Apt unlike sign pair - else - ptTrig = 0x200;// mask for Lpt unlike sign pair + ptTrig = 0x04;// mask for Lpt unlike sign pair if (esd->GetTriggerMask() & ptTrig) NbTrigger++; @@ -463,14 +496,16 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven hPtResonance->Fill(fVtot.Pt()); // match with trigger - if (muonTrack->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig)) EventInMassMatch++; + if (muonTrack2->GetMatchTrigger()>=0 && (esd->GetTriggerMask() & ptTrig)) EventInMassMatch++; } } // 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); @@ -562,7 +597,17 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven cout << "Chi2Cut for muon tracks = " << Chi2Cut << endl; cout << "PtCutMin for muon tracks = " << PtCutMin << endl; cout << "PtCutMax for muon tracks = " << PtCutMax << endl; - cout << "Entries (unlike sign dimuons) in the mass range ["<GetEntries(); + + if (hInvMassAll->GetEntries() > 0) { + hInvMassAll->Fit("gaus","q0"); + TF1* f1 = hInvMassAll->GetFunction("gaus"); + cout << Form(". Rough sigma = %7.2f MeV/c2",f1->GetParameter(2)*1000.0); + } + + cout << endl << "Entries (unlike sign dimuons) in the mass range ["<