smaller test macro for ESD track matching
authorjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Aug 2008 17:05:00 +0000 (17:05 +0000)
committerjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Aug 2008 17:05:00 +0000 (17:05 +0000)
EMCAL/macros/TestESD.C [new file with mode: 0644]

diff --git a/EMCAL/macros/TestESD.C b/EMCAL/macros/TestESD.C
new file mode 100644 (file)
index 0000000..4e6fe2c
--- /dev/null
@@ -0,0 +1,108 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+//Root include files 
+#include <Riostream.h>
+#include <TFile.h>
+#include <TChain.h>
+#include <TParticle.h>
+#include <TNtuple.h>
+#include <TCanvas.h>
+#include <TObjArray.h>
+#include <TSystem.h>
+#include <TString.h>
+#include <TH1F.h>
+#include <TVector.h>
+#include <TParticle.h>
+#include <TRefArray.h>
+#include <TArrayS.h>
+
+//AliRoot include files 
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliESDCaloCluster.h"
+#include "AliESDCaloCells.h"
+#include "AliPID.h"
+#include "AliLog.h"
+
+#endif
+
+void TestESD() {
+
+  TFile* f = new TFile("AliESDs.root");
+  TTree* esdTree = (TTree*)f->Get("esdTree");
+  
+  AliESDEvent* esd = new AliESDEvent();
+  esd->ReadFromTree(esdTree);
+
+  Int_t nEvt = esdTree->GetEntries();
+  Float_t pos[3] = {0.,0.,0.};
+
+  for(Int_t iev = 0; iev < nEvt; iev++) {
+    cout << "Event: " << iev+1 << "/" << nEvt << endl;
+    esdTree->GetEvent(iev);
+
+    TRefArray* caloClusters = new TRefArray();
+    esd->GetEMCALClusters(caloClusters);
+
+    //get reconstructed vertex position
+    Double_t vertex_position[3];
+    esd->GetVertex()->GetXYZ(vertex_position);
+
+    //loop over clusters
+    Int_t nclus = caloClusters->GetEntries();
+    for (Int_t icl = 0; icl < nclus; icl++) {
+
+      AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
+      Float_t energy = clus->E();
+      clus->GetPosition(pos);
+      TVector3 vpos(pos[0],pos[1],pos[2]);
+      TLorentzVector p;
+      clus->GetMomentum(p,vertex_position);
+      Double_t cphi = vpos.Phi();
+      Double_t ceta = vpos.Eta();
+
+      Int_t nMatched = clus->GetNTracksMatched();
+      Int_t trackIndex = clus->GetTrackMatched();
+      Int_t nLabels = clus->GetNLabels();
+      Int_t labelIndex = clus->GetLabel();
+
+      Int_t nCells = clus->GetNCells();
+
+      //For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME
+
+      cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: " 
+          << cphi << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched 
+          << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: " 
+          << labelIndex << endl;
+
+      if(trackIndex >= 0) {
+       AliESDtrack* track = esd->GetTrack(trackIndex);
+       Double_t tphi = track->GetOuterParam()->Phi();
+       Double_t teta = track->GetOuterParam()->Eta();
+       Double_t tmom = track->GetOuterParam()->P();
+       cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
+
+       Double_t deta = teta - ceta;
+       Double_t dphi = tphi - cphi;
+       if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
+       if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
+       Double_t dR = sqrt(dphi*dphi + deta*deta);
+
+       Double_t pOverE = tmom/energy;
+
+       if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
+         cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
+       }
+
+      }
+
+    }
+
+
+  }
+
+
+
+}