track matching macros from Alberto
authorjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Oct 2006 05:44:51 +0000 (05:44 +0000)
committerjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Oct 2006 05:44:51 +0000 (05:44 +0000)
EMCAL/macros/trackMatching/DrawComparison.C [new file with mode: 0755]
EMCAL/macros/trackMatching/FindMatches.C [new file with mode: 0755]
EMCAL/macros/trackMatching/MatchComparison.C [new file with mode: 0755]
EMCAL/macros/trackMatching/SaveTrueMatchesSimple.C [new file with mode: 0755]

diff --git a/EMCAL/macros/trackMatching/DrawComparison.C b/EMCAL/macros/trackMatching/DrawComparison.C
new file mode 100755 (executable)
index 0000000..5b6db4e
--- /dev/null
@@ -0,0 +1,42 @@
+//
+// This macro reads comparison files and makes plots.
+//
+
+void DrawComparison
+(const char *fileName = "match-comparison.root")
+{
+       TH1D *hg = new TH1D("hg", "Good matches / true matches", 20, 0.0,  10.0);
+       TH1D *hf = new TH1D("hf", "Fake matches / found matches", 20, 0.0,  10.0);
+       
+       // 
+       // Open file
+       //
+       TFile *file = TFile::Open(fileName);
+       if (!file) return;
+       
+       TH1D *hgood = (TH1D*)file->Get("hgood");
+       TH1D *hfake = (TH1D*)file->Get("hfake");
+       TH1D *htrue = (TH1D*)file->Get("htrue");
+       TH1D *hfound = (TH1D*)file->Get("hfound");
+
+       gROOT->SetStyle("Plain");
+       gStyle->SetOptStat(0);
+       
+       hg->Divide(hgood, htrue, 100.0, 1.0, "b");
+       hf->Divide(hfake, hfound, 100.0, 1.0, "b");
+       
+       TCanvas *c = new TCanvas("c", "", 0, 0, 800, 600);
+       
+       hg->SetMarkerStyle(21);
+       hf->SetMarkerStyle(25);
+       
+       hg->SetXTitle("p_{T} (GeV/c)");
+
+       hg->GetXaxis()->SetRangeUser(0.0, 6.0);
+       hf->GetXaxis()->SetRangeUser(0.0, 6.0);
+       
+       hg->SetMaximum(120.0);
+       hg->SetMinimum(0.0);
+       hg->Draw("PE1");
+       hf->Draw("PE1same");
+}      
diff --git a/EMCAL/macros/trackMatching/FindMatches.C b/EMCAL/macros/trackMatching/FindMatches.C
new file mode 100755 (executable)
index 0000000..54b17e5
--- /dev/null
@@ -0,0 +1,105 @@
+//
+// This macro performs the matching operation between ESD tracks 
+// and EMCAL clusters (stored as AliESDCaloCluster objects).
+// Due to the necessity to know the magnetic field, it is supposed that
+// this macro runs in the directory where event files are stored
+// (galice.root --> recovery magnetic field), and saves its output in
+// a directory specified in the argument (default = work directory)
+// 
+
+void FindMatches(const char *fileOut = "matchESD.root")
+{
+       //
+       // Initialize AliRun manager.
+       //
+       if (gAlice) {
+               delete gAlice;
+               gAlice = 0;
+       }
+       
+       //
+       // Initialize AliRunLoader, in order
+       // to read magnetic field from simulation.
+       //
+       AliRunLoader *rl = AliRunLoader::Open("galice.root");
+       if (!rl) return;
+       rl->LoadgAlice();
+       gAlice = rl->GetAliRun();
+       AliMagF *magf = gAlice->Field();
+       Bool_t constField = (magf->Type() == 1);
+       AliTracker::SetFieldMap(magf, constField);
+       
+       //
+       // Open ESD file and recoveries TTree of ESD objects.
+       //
+       TFile *esdFile = new TFile("AliESDs.root");
+       TTree *esdTree = (TTree*)esdFile->Get("esdTree");
+       AliESD *esd = 0;
+       esdTree->SetBranchAddress("ESD", &esd);
+       Long64_t nEvents = esdTree->GetEntries();
+       
+       //
+       // Set this important flag to the tracker.
+       // This example works with and AliESD file already saved
+       // after the tracking procedure done during AliReconstruction::Run().
+       // All tracks saved in such a file, are already propagated to the vertex,
+       // while EMCAL matching needs to be done using track status in the outermost
+       // available position.
+       // Then, we must use the "outer" parameters in the ESD track, and we must
+       // tell to the track to copy the outer parameters from its AliESDtrack seed.
+       // 
+       AliEMCALTrack::SetUseOuterParams(kTRUE);
+       
+       //
+       // Instantiate match finder, and set some cuts.
+       // Distances are in centimeters, angles in degrees.
+       //
+       AliEMCALTracker *mf = new AliEMCALTracker;
+       mf->SetCutX(5.0);
+       mf->SetCutY(5.0);
+       mf->SetCutZ(5.0);
+       mf->SetCutAlpha(-50., 50.);  // --> i.e. exclude this cut
+       mf->SetCutAngle(10.);
+       mf->SetMaxDistance(5.0);
+       
+       //
+       // Define how to manage energy loss correction.
+       // Actually, these corrections are exlcluded.
+       //
+       // mf->SetCorrection(0.0604557, 34.5437); // at the moment, we exclude energy loss correction
+       mf->SetTrackCorrectionMode("NONE");
+       mf->SetNumberOfSteps(0);
+       //TGeoManager::Import("misaligned_geometry.root");      
+       
+       //
+       // Before starting looping on files, the output file is created.
+       // It will be structurally identical to the ESD source tree, with the 
+       // only difference that here the "fEMCALindex" datamember of the ESD tracks
+       // will have been set to a meaningful value.
+       //
+       TFile *outFile = TFile::Open(fileOut, "RECREATE");
+       TTree *outTree = new TTree("esdTree", "ESD with matched clusters");
+       outTree->Branch("ESD", "AliESD", &esd);
+       
+       //
+       // Loop on events.
+       // The whole track matching procedure is compiled in the 
+       // method "PropagateBack" which accepts and input ESD collection.
+       // This method modifies the passed object then, the same object is linked
+       // to source tree and target tree branch.
+       //
+       Int_t nTracks, nStep1, nStep2, nSaved;
+       for (Long64_t iev = 0; iev < nEvents; iev++) {
+               cout << "Finding matches in event " << iev + 1 << "/" << nEvents << endl;
+               esdTree->GetEntry(iev);
+               mf->PropagateBack(esd);
+               outTree->Fill();
+       }
+       
+       // 
+       // Save processed tree and close output file.
+       //
+       outFile->cd();
+       outTree->Write("esdTree", TObject::kOverwrite);
+       outFile->Close();
+}
diff --git a/EMCAL/macros/trackMatching/MatchComparison.C b/EMCAL/macros/trackMatching/MatchComparison.C
new file mode 100755 (executable)
index 0000000..64fc2db
--- /dev/null
@@ -0,0 +1,132 @@
+//
+// This macros compares the matches found by algorithm with the true matches
+// found with the "SaveTrueMatchesSimple.C" macro.
+// It saves 4 histogram, which contain the Pt distribution of:
+//   - all found matches
+//   - all correctly found matches
+//   - all wrong (fake) found matches
+//   - all true matches
+// which will then be availabel for computing efficiency and contamination.
+//
+
+class match_t
+{
+       public:
+       
+       Int_t    label;       // GEANT label of particle
+       Int_t    indexT;      // index of track in ESD collection
+       Int_t    indexC;      // index of cluster in ESD collection
+       Double_t p[3];        // track momentum
+       Double_t v[3];        // track vertex
+};
+
+void MatchComparison()
+{
+       //
+       // Initialize AliRun manager
+       //
+       if (gAlice) {
+               delete gAlice;
+               gAlice = 0;
+       }
+       
+       //
+       // Initialize run loader and load Kinematics
+       //
+       AliRunLoader *runLoader = AliRunLoader::Open("galice.root");
+       if (!runLoader) return;
+       runLoader->LoadgAlice();
+       gAlice = runLoader->GetAliRun();
+       runLoader->LoadKinematics();
+       
+       //
+       // Initialize histograms with their error computation
+       //
+       TH1D *hgood = new TH1D("hgood", "Well matched tracks", 20, 0.0,  10.0);
+       TH1D *hfake = new TH1D("hfake", "Fake matched tracks", 20, 0.0,  10.0);
+       TH1D *htrue = new TH1D("htrue", "True matches", 20, 0.0,  10.0);
+       TH1D *hfound = new TH1D("hfound", "Found matches", 20, 0.0,  10.0);
+       hgood->Sumw2();
+       hfake->Sumw2();
+       htrue->Sumw2();
+       hfound->Sumw2();
+       
+       //
+       // Open file containing true matches,
+       // retrieve the Tree and link to a cursor.
+       //
+       TFile *fileTrue = TFile::Open("true-matches.root");
+       match_t trueMatch;
+       
+       //
+       // Open file of found matches,
+       // link the modified ESD container.
+       //
+       TFile *fileFound = TFile::Open("matchESD.root");
+       TTree *treeFound = (TTree*)fileFound->Get("esdTree");
+       AliESD *esd = 0;
+       treeFound->SetBranchAddress("ESD", &esd);
+       Long64_t nEvents = treeFound->GetEntries();
+       
+       //
+       // Loop on all events
+       //
+       Int_t im, it, ic, nTrueMatches, nTracks;
+       Int_t label, trkLabel, cluLabel;
+       for (Long64_t iev = 0; iev < nEvents; iev++) {
+               
+               // get true matches tree of given event
+               TTree *treeTrue = (TTree*)fileTrue->Get(Form("tm_%d", iev));
+               treeTrue->SetBranchAddress("matches", &trueMatch);
+               nTrueMatches = treeTrue->GetEntries();
+               
+               // set TTree pointers to selected event
+               runLoader->GetEvent(iev);
+               treeFound->GetEntry(iev);
+               AliStack *stack = runLoader->Stack();
+               nTracks = esd->GetNumberOfTracks();
+               
+               // read all true pairs
+               for (im = 0; im < nTrueMatches; im++) {
+                       treeTrue->GetEntry(im);
+                       AliESDtrack *track = esd->GetTrack(trueMatch.indexT);
+                       label = TMath::Abs(track->GetLabel());
+                       TParticle *p = stack->Particle(label);
+                       htrue->Fill(p->Pt());
+               }
+               
+               // compare found matches
+               for (Int_t it = 0; it < nTracks; it++) {
+                       AliESDtrack *track = esd->GetTrack(it);
+                       ic = track->GetEMCALcluster();
+                       if (ic == -999999999) continue;
+                       ic = TMath::Abs(ic);
+                       AliESDCaloCluster *cl = esd->GetCaloCluster(ic);
+                       if (!cl) continue;
+                       trkLabel = track->GetLabel();
+                       cluLabel = cl->GetPrimaryIndex();
+                       if (trkLabel == cluLabel && trkLabel > 0) {
+                               TParticle *p = stack->Particle(TMath::Abs(trkLabel));
+                               hgood->Fill(p->Pt());
+                               hfound->Fill(p->Pt());
+                       }
+                       else  {
+                               TParticle *p = stack->Particle(TMath::Abs(trkLabel));
+                               hfake->Fill(p->Pt());
+                               hfound->Fill(p->Pt());
+                       }
+               }
+       }
+       
+       cout << "True matches : " << htrue->GetEntries() << endl;
+       cout << "Found matches: " << hfound->GetEntries() << endl;
+       cout << "Good matches : " << hgood->GetEntries() << endl;
+       cout << "Fake matches : " << hfake->GetEntries() << endl;
+       
+       TFile *fout = TFile::Open("match-comparison.root", "RECREATE");
+       hgood->Write();
+       hfake->Write();
+       htrue->Write();
+       hfound->Write();
+       fout->Close();
+}
diff --git a/EMCAL/macros/trackMatching/SaveTrueMatchesSimple.C b/EMCAL/macros/trackMatching/SaveTrueMatchesSimple.C
new file mode 100755 (executable)
index 0000000..fec167e
--- /dev/null
@@ -0,0 +1,128 @@
+//
+// This macro generates a simple TTree containing
+// all true matches collected from one event.
+//
+// For each match it is stored:
+//  - label ID of the track
+//  - label ID of the EMCAL cluster
+//  - momentum (3D) of the track
+//  - vertex (3D) of the track
+//
+// All tracks which come from kinks are excluded 
+// from this computation, and fake tracks are collected.
+// If desired, it is possible to change these settings 
+// operating on the Bool_t variables listed below.
+//  
+
+Bool_t rejectKinks = kTRUE;     // switch to TRUE to include kinks in computation
+Bool_t rejectFakes = kFALSE;    // switch to TRUE to include fake tracks in computation
+
+class match_t
+{
+       public:
+       
+       Int_t    label;       // GEANT label of particle
+       Int_t    indexT;      // index of track in ESD collection
+       Int_t    indexC;      // index of cluster in ESD collection
+       Double_t p[3];        // track momentum
+       Double_t v[3];        // track vertex
+};
+
+//
+// Read AliESDs.root file and saves all true pairs of trak-cluster.
+//
+void SaveTrueMatchesSimple_compiled(const char *outFileName)
+{
+       //
+       // open ESD file, retrieve tree and link to branch cursor
+       //
+       TFile *srcFile = TFile::Open("AliESDs.root");
+       if (!srcFile) return;
+       TTree *srcTree = (TTree*)srcFile->Get("esdTree");
+       if (!srcTree) return;
+       AliESD *esd = 0;
+       srcTree->SetBranchAddress("ESD", &esd);
+       Long64_t nEvents = srcTree->GetEntries();
+       
+       //
+       // Open output file and create output tree
+       //
+       TFile *outFile = new TFile(outFileName, "RECREATE");
+       match_t match;
+
+       //
+       // Loop on events and store true matches
+       //
+       Bool_t isKink;
+       Int_t label, count, nTracks, firstCluster, lastCluster;
+       for (Long64_t iev = 0; iev < nEvents; iev++) {
+
+               srcTree->GetEntry(iev);
+               cout << "Event " << iev + 1 << " of " << nEvents << ": " << endl;
+               
+               nTracks = esd->GetNumberOfTracks();
+               firstCluster = esd->GetFirstEMCALCluster();
+               lastCluster = esd->GetFirstEMCALCluster() + esd->GetNumberOfEMCALClusters();
+               cout << "Tracks found      : " << nTracks << endl;
+               cout << "EMC clusters found: " << lastCluster - firstCluster << endl;
+               
+               // create new matches tree
+               TTree *outTree = new TTree(Form("tm_%d", iev), Form("True matches from event %d", iev));
+               outTree->Branch("matches", &match, "label/I:indexT/I:indexC/I:p[3]/D:v[3]/D");
+               
+               // external loop on tracks
+               for (Int_t it = 0; it < nTracks; it++) {
+                       AliESDtrack *esdTrack = esd->GetTrack(it);
+                       // start check to reject kinks
+                       if (rejectKinks) {
+                               isKink = kFALSE;
+                               for (Int_t i = 0; i < 3; i++) {
+                                       if (esdTrack->GetKinkIndex(i) > 0) isKink = kTRUE;
+                               }
+                               if (isKink) continue;
+                       }
+                       // get track GEANT label (to be checked for match)
+                       label = esdTrack->GetLabel();
+                       if (rejectFakes) {
+                               if (label < 0) continue;
+                       }
+                       else {
+                               label = TMath::Abs(label);
+                       }
+                       // store track data into candidate match variable
+                       // anyway, it will be stored in the tree only
+                       // if it matches a cluster
+                       match.indexT = it;
+                       esdTrack->GetPxPyPz(match.p);
+                       esdTrack->GetXYZ(match.v);
+                       // internal loop on clusters
+                       // a counter counts how many true matches are
+                       // generated for the same track
+                       count = 0;
+                       for (Int_t ic = firstCluster; ic < lastCluster; ic++) {
+                               AliESDCaloCluster *cl = esd->GetCaloCluster(ic);
+                               // reject pseudo-clusters & unmatched clusters
+                               if (cl->GetClusterType() != AliESDCaloCluster::kClusterv1) continue;
+                               if (cl->GetPrimaryIndex() != label) continue;
+                               // if the method reaches this point, we
+                               // have found a match to be stored
+                               match.label = label;
+                               match.indexC = ic;
+                               outTree->Fill();
+                               count++;
+                       }
+                       // alert for multiple matches
+                       if (count > 0) {
+                               cout << "Found " << count << " clusters which match track " << it << " in ESD";
+                               if (count > 1) cout << " --> MULTIPLE MATCH";
+                               cout << endl;
+                       }
+               } // end loop on tracks
+               
+               outFile->cd();
+               outTree->Write();
+               delete outTree;
+       }
+       
+       outFile->Close();
+}