ITS/testITSMultReco.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 #include <Riostream.h>
   4 #include <TFile.h>
   5 #include <TTree.h>
   6 #include <TBranch.h>
   7 #include <TGeoManager.h>
   8
   9 #include "AliRunLoader.h"
  10 #include "AliESD.h"
  11 #include "AliRun.h"
  12
  13 #include "AliITSgeom.h"
  14 #include "AliITSLoader.h"
  15 #include "AliITSMultReconstructor.h"
  16
  17 #endif
  18
  19 void testITSMultReco(Char_t* dir = ".") {
  20
  21   Char_t str[256];
  22
  23   // ########################################################
  24   // defining pointers
  25   AliRunLoader* runLoader;
  26   TFile* esdFile = 0;
  27   TTree* esdTree = 0;
  28   AliESD* esd = 0;
  29
  30   // #########################################################
  31   // get runloader
  32
  33   if (gAlice) {
  34     delete gAlice->GetRunLoader();
  35     delete gAlice;
  36     gAlice=0;
  37   }
  38
  39   sprintf(str,"%s/galice.root",dir);
  40   runLoader = AliRunLoader::Open(str);
  41   if (runLoader == 0x0) {
  42     cout << "Can not open session"<<endl;
  43     return;
  44   }
  45   // get geometry (here geometry.root is used, change it if needed)
  46   if (!gGeoManager) {
  47     sprintf(str,"%s/geometry.root",dir);
  48     TGeoManager::Import(str);
  49     if(!gGeoManager) {
  50       cout << "Can not access the geometry file"<<endl;
  51       return;
  52     }
  53   }
  54
  55
  56   // #########################################################
  57   // open esd file and get the tree
  58
  59   // close it first to avoid memory leak
  60   if (esdFile)
  61     if (esdFile->IsOpen())
  62       esdFile->Close();
  63
  64   sprintf(str,"%s/AliESDs.root",dir);
  65   esdFile = TFile::Open(str);
  66   esdTree = (TTree*)esdFile->Get("esdTree");
  67   TBranch * esdBranch = esdTree->GetBranch("ESD");
  68   esdBranch->SetAddress(&esd);
  69
  70
  71   // #########################################################
  72   // setup its stuff
  73
  74   AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader");
  75   if (!itsLoader) {
  76     cout << " Can't get the ITS loader!" << endl;
  77     return ;
  78   }
  79   AliITSgeom* itsGeo=itsLoader->GetITSgeom();
  80   itsLoader->LoadRecPoints("read");
  81
  82   // #########################################################
  83   AliITSMultReconstructor* multReco = new AliITSMultReconstructor();
  84   multReco->SetGeometry(itsGeo);
  85
  86   // #########################################################
  87   // getting number of events
  88
  89   Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
  90   Int_t nESDEvents = esdBranch->GetEntries();
  91
  92   if (nEvents!=nESDEvents) {
  93     cout << " Different number of events from runloader and esdtree!!!"
  94          << nEvents << " / " << nESDEvents << endl;
  95     return;
  96   }
  97
  98   // ########################################################
  99   // loop over number of events
 100   cout << nEvents << " event(s) found in the file set" << endl;
 101   for(Int_t i=0; i<nEvents; i++) {
 102
 103     cout << "-------------------------" << endl << " event# " << i << endl;
 104
 105     runLoader->GetEvent(i);
 106     esdBranch->GetEntry(i);
 107
 108     // ########################################################
 109     // get the EDS vertex
 110     const AliESDVertex* vtxESD = esd->GetVertex();
 111     Double_t vtx[3];
 112     vtxESD->GetXYZ(vtx);
 113     Float_t esdVtx[3];
 114     esdVtx[0] = vtx[0];
 115     esdVtx[1] = vtx[1];
 116     esdVtx[2] = vtx[2];
 117
 118     ///#########################################################
 119     // get ITS clusters
 120     TTree* itsClusterTree = itsLoader->TreeR();
 121     if (!itsClusterTree) {
 122       cerr<< " Can't get the ITS cluster tree !\n";
 123       return;
 124     }
 125     multReco->SetHistOn(kTRUE);
 126     multReco->Reconstruct(itsClusterTree, esdVtx, esdVtx);
 127
 128     cout <<" >>>> Number of tracklets: "<<multReco->GetNTracklets()<<endl;
 129     for (Int_t t=0; t<multReco->GetNTracklets(); t++) {
 130
 131       cout << "  tracklet " << t
 132            << " , theta = " << multReco->GetTracklet(t)[0]
 133            << " , phi = " << multReco->GetTracklet(t)[1]
 134            << " , DeltaPhi = " << multReco->GetTracklet(t)[2]<< endl;
 135     }
 136     cout <<" >>>> Number of single layer 1 clusters: "<<multReco->GetNSingleClusters()<<endl;
 137     for (Int_t t=0; t<multReco->GetNSingleClusters(); t++) {
 138
 139       cout << "  cluster " << t
 140            << " , theta = " << multReco->GetCluster(t)[0]
 141            << " , phi = " << multReco->GetCluster(t)[1] << endl;
 142     }
 143
 144   }
 145
 146   TFile* fout = new TFile("out.root","RECREATE");
 147
 148   multReco->SaveHists();
 149   fout->Write();
 150   fout->Close();
 151
 152
 153 }