/* $Id$ */
-// A. De Falco, H. Woehri, INFN Cagliari, July 2006
-// This macro reads the generation/reconstruction files in the
-// input directory (dirname),
-// builds a tree that contains, for each event, an array of muons and dimuons
-// (TClonesArrays of AliMUONTrackLight and AliMUONPairLight objects)
-// and writes the tree in an output file (outFileName)
-// Note that if the path for the output file is not explicitly specified,
-// it will be written in the directory containing the generation/reconstruction
+/// \ingroup macros
+/// \file DecodeRecoCocktail.C
+/// \brief add brief description
+///
+/// \author A. De Falco, H. Woehri, INFN Cagliari, July 2006
+///
+/// This macro reads the generation/reconstruction files in the
+/// input directories (recodir and simdir),
+/// builds a tree that contains, for each event, an array of muons and dimuons
+/// (TClonesArrays of AliMUONTrackLight and AliMUONPairLight objects)
+/// and writes the tree in an output file (outFileName)
+/// Note that if the path for the output file is not explicitly specified,
+/// it will be written in the directory containing the generation/reconstruction
+/// 27-Nov-2006: modified by in order to loop on files
+///
+/// 13 Nov 2007:
+/// Updated this macro to work with new version of AliMUONRecoCheck. Also, we are
+/// now fetching reconstructed track data from ESD and not the track tree, because
+/// the AliMUONTrack objects for reconstructed tracks are no longer stored on disk.
+/// - Artur Szostak <artursz@iafrica.com>
-#include <iostream>
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>
#include <TLorentzVector.h>
#include <TParticle.h>
#include <TSystem.h>
-#include <AliMUONRecoCheck.h>
-#include <AliMUONTrack.h>
-#include <AliMUONTrackParam.h>
-#include <AliRunLoader.h>
+#include "AliMUONRecoCheck.h"
+#include "AliMUONTrack.h"
#include "AliMUONTrackLight.h"
#include "AliMUONPairLight.h"
+#include "AliMUONVTrackStore.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliESDMuonTrack.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+/*TODO: need to update this with changes made to ITS
+#include "AliITSVertexerPPZ.h"
+#include "AliITSLoader.h"
+*/
+#endif
-void DecodeRecoCocktail(char* dirname=".", char* outFileName = "MuonLight.root"){
- const char *startingDir = gSystem->pwd();
- gSystem->cd(dirname);
- TFile *fout = new TFile(outFileName,"recreate");
- AliMUONRecoCheck *rc = new AliMUONRecoCheck("galice.root");
- AliRunLoader *runLoader = rc->GetRunLoader();
-
+void DecodeRecoCocktail(
+ char* recodir=".", // The directory containing galice.root for reconstructed data.
+ char* simdir="generated/", // The directory containing galice.root for simulated data.
+ char* outFileName = "MuonLight.root" // The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
+ )
+{
+/// \param recodir The directory containing galice.root for reconstructed data.
+/// \param simdir The directory containing galice.root for simulated data.
+/// \param outFileName The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
+
+ char startingDir[200];
+ sprintf (startingDir,"%s",gSystem->pwd());
+ gSystem->cd(recodir);
+
TClonesArray *muonArray = new TClonesArray("AliMUONTrackLight",10);
TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
TTree *treeOut = new TTree("tree","tree");
+ TFile *fout = new TFile(outFileName,"recreate");
fout->cd();
+
treeOut->Branch("muons",&muonArray);
treeOut->Branch("dimuons",&dimuonArray);
- runLoader->LoadKinematics("READ");;
- Int_t nev = runLoader->GetNumberOfEvents();
+ AliMUONRecoCheck *rc = new AliMUONRecoCheck("AliESDs.root", simdir);
+
+ /*TODO: need to update this with changes made to ITS
+ AliITSLoader* ITSloader = (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
+ AliITSVertexerPPZ *dovert = 0;
+ if (ITSloader) {
+ dovert = new AliITSVertexerPPZ("default",0,0);
+ // dovert->SetDebug(0);
+ dovert->SetDiffPhiMax(0.05);
+ dovert->SetWindow(3.);
+ }
+ AliESDVertex *vert = 0;
+ */
TLorentzVector v;
+
+ Int_t nev = rc->NumberOfEvents();
for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
- printf ("Event %d of %d\n",ievent,nev);
+
+ /*TODO: need to update this with changes made to ITS
+ runLoaderSim->GetHeader();
+ if (ITSloader) {
+ vert = dovert->FindVertexForCurrentEvent(ievent);
+ }
+ */
+
+ //printf ("Event %d of %d\n",ievent,nev);
muonArray->Clear(); // clean muon and dimuon arrays
dimuonArray->Clear();
- runLoader->GetEvent(ievent);
- rc->ResetTracks();
- rc->MakeTrackRef(); // make reconstructable tracks
- TClonesArray * trackRecoArray = rc->GetTrackReco();
- TClonesArray * trackRefArray = rc->GetMuonTrackRef();
- Int_t nTrackReco = trackRecoArray->GetEntriesFast();
- Int_t nreftracks = 0;
- for (Int_t itrRec = 0; itrRec<nTrackReco; itrRec++) {
- // assign parameters concerning the reconstructed tracks
- AliMUONTrack *trackReco = (AliMUONTrack *)trackRecoArray->At(itrRec);
- AliMUONTrackLight *muLight = new AliMUONTrackLight();
- AliMUONTrackParam *trPar = trackReco->GetTrackParamAtVertex();
- muLight->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
- muLight->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
- muLight->SetTriggered(trackReco->GetMatchTrigger());
- Double_t xyz[3] = { trPar->GetNonBendingCoor(),
- trPar->GetBendingCoor(),
- trPar->GetZ()};
- muLight->SetVertex(xyz);
- // find the reference track and store further information
- TParticle *part = muLight->FindRefTrack(trackReco,trackRefArray,runLoader);
- if (part) {
+ AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(ievent, kFALSE); // Use tracks from actual reconstruction.
+ //AliMUONVTrackStore* trackRefs = rc->ReconstructibleTracks(ievent); // Only use reconstructible reference tracks.
+ AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent);
+ AliStack* pstack = (const_cast<AliMCEventHandler*>(rc->GetMCEventHandler()))->MCEvent()->Stack();
+
+ // loop over ESD tracks
+ Int_t nreftracks = 0;
+ const AliESDEvent* esd = rc->GetESDEvent();
+ Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
+ for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+
+ AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
+
+ // skip ghosts
+ if (!esdTrack->ContainTrackerData()) continue;
+
+ // find the corresponding MUON track
+ AliMUONTrack* trackReco = (AliMUONTrack*) recoTracks->FindObject(esdTrack->GetUniqueID());
+
+ // try to match the reconstructed track with a simulated one
+ Int_t nMatchClusters = 0;
+ AliMUONTrack* matchedTrackRef = rc->FindCompatibleTrack(*trackReco, *trackRefs, nMatchClusters, kFALSE, 10.);
+
+ if (matchedTrackRef) {
+
+ //store new referenced track in the muonArray
+ AliMUONTrackLight* muLight = new ((*muonArray)[nreftracks++]) AliMUONTrackLight();
+
+ // assign parameters concerning the reconstructed tracks
+ muLight->FillFromESD(esdTrack);
+ // muLight->FillFromAliMUONTrack(trackReco);
+
+ // store further information related to the simulated track
+ muLight->SetTrackPythiaLine(matchedTrackRef->GetUniqueID());
+ TParticle *part = pstack->Particle(matchedTrackRef->GetUniqueID());
+ muLight->SetTrackPDGCode(part->GetPdgCode());
v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
muLight->SetPGen(v);
- muLight->FillMuonHistory(runLoader, part);
- muLight->PrintInfo("A");
- //store the referenced track in the muonArray:
- TClonesArray &muons = *muonArray;
- new (muons[nreftracks++]) AliMUONTrackLight(*muLight);
+ muLight->FillMuonHistory(pstack, part);
+ // muLight.PrintInfo("A");
+
}
- }
+
+ } // end esd track
// now loop over muon pairs to build dimuons
Int_t nmuons = muonArray->GetEntriesFast();
AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
- AliMUONPairLight *dimuLight = new AliMUONPairLight();
- dimuLight->SetMuons(*mu1, *mu2);
- // if(dimuLight->GetCreationProcess() == 2){
- // printf ("#dimuon = %d (%d, %d) \n", ndimuons, itrRec1, itrRec2);
- // dimuLight->PrintInfo("A");
- // }
- //store the referenced track in the dimuonArray:
+ AliMUONPairLight dimuLight;
+ dimuLight.SetMuons(*mu1, *mu2);
TClonesArray &dimuons = *dimuonArray;
- new (dimuons[ndimuons++]) AliMUONPairLight(*dimuLight);
+ new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
}
}
- Int_t ndimu2 = dimuonArray->GetEntriesFast();
- // printf ("dimuonArray has %d entries\n",ndimu2);
- // dimuonArray->Dump();
- // printf ("filling tree\n");
- // // fill the tree
- treeOut->Fill();
- // printf ("done\n");
- }
+ treeOut->Fill();
+ }
+
fout->cd();
treeOut->Write();
gSystem->cd(startingDir);
+ fout->Close();
+ delete muonArray;
+ delete dimuonArray;
+ delete treeOut;
+ delete fout;
+ delete rc;
}
+