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
   8 #include "AliRunLoader.h"
   9 #include "AliESD.h"
  10 #include "AliRun.h"
  11
  12 #include "AliITS.h"
  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   // setup galice and 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   runLoader->LoadgAlice();
  46
  47   gAlice = runLoader->GetAliRun();
  48   runLoader->LoadKinematics();
  49   runLoader->LoadHeader();
  50
  51   // #########################################################
  52   // open esd file and get the tree
  53
  54   // close it first to avoid memory leak
  55   if (esdFile)
  56     if (esdFile->IsOpen())
  57       esdFile->Close();
  58
  59   sprintf(str,"%s/AliESDs.root",dir);
  60   esdFile = TFile::Open(str);
  61   esdTree = (TTree*)esdFile->Get("esdTree");
  62   TBranch * esdBranch = esdTree->GetBranch("ESD");
  63   esdBranch->SetAddress(&esd);
  64
  65
  66   // #########################################################
  67   // setup its stuff
  68
  69   AliITS* its=(AliITS*)runLoader->GetAliRun()->GetDetector("ITS");
  70   if (!its) {
  71     cout << " Can't get the ITS!" << endl;
  72     return ;
  73   }
  74   AliITSgeom* itsGeo=its->GetITSgeom();
  75   if (!itsGeo) {
  76     cout << " Can't get the ITS geometry!" << endl;
  77     return ;
  78   }
  79   AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader");
  80   if (!itsLoader) {
  81     cout << " Can't get the ITS loader!" << endl;
  82     return ;
  83   }
  84   itsLoader->LoadRecPoints("read");
  85
  86   // #########################################################
  87   AliITSMultReconstructor* multReco = new AliITSMultReconstructor();
  88   multReco->SetGeometry(itsGeo);
  89
  90   // #########################################################
  91   // getting number of events
  92
  93   Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
  94   Int_t nESDEvents = esdBranch->GetEntries();
  95
  96   if (nEvents!=nESDEvents) {
  97     cout << " Different number of events from runloader and esdtree!!!"
  98          << nEvents << " / " << nESDEvents << endl;
  99     return;
 100   }
 101
 102   // ########################################################
 103   // loop over number of events
 104   cout << nEvents << " event(s) found in the file set" << endl;
 105   for(Int_t i=0; i<nEvents; i++) {
 106
 107     cout << "-------------------------" << endl << " event# " << i << endl;
 108
 109     runLoader->GetEvent(i);
 110     esdBranch->GetEntry(i);
 111
 112     // ########################################################
 113     // get the EDS vertex
 114     const AliESDVertex* vtxESD = esd->GetVertex();
 115     Double_t vtx[3];
 116     vtxESD->GetXYZ(vtx);
 117     Float_t esdVtx[3];
 118     esdVtx[0] = vtx[0];
 119     esdVtx[1] = vtx[1];
 120     esdVtx[2] = vtx[2];
 121
 122     ///#########################################################
 123     // get ITS clusters
 124     TTree* itsClusterTree = itsLoader->TreeR();
 125     if (!itsClusterTree) {
 126       cerr<< " Can't get the ITS cluster tree !\n";
 127       return;
 128     }
 129     multReco->SetHistOn(kTRUE);
 130     multReco->Reconstruct(itsClusterTree, esdVtx, esdVtx);
 131
 132
 133     for (Int_t t=0; t<multReco->GetNTracklets(); t++) {
 134
 135       cout << "  tracklet " << t
 136            << " , theta = " << multReco->GetTracklet(t)[0]
 137            << " , phi = " << multReco->GetTracklet(t)[1] << endl;
 138     }
 139
 140   }
 141
 142   TFile* fout = new TFile("out.root","RECREATE");
 143
 144   multReco->SaveHists();
 145   fout->Write();
 146   fout->Close();
 147
 148
 149 }