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