Tracking V1 ported to the HEAD
[u/mrichter/AliRoot.git] / ITS / AliITSStoreFindableTracksCompiled.C
1 #include <fstream.h>
2
3 #include <TClassTable.h>
4 #include <TFile.h>
5 #include <TTree.h>
6 #include <TString.h>
7 #include <TClonesArray.h>
8 #include <TParticle.h>
9
10 #include "AliRun.h"
11 #include "AliITS.h"
12 #include "AliITSgeom.h"
13 #include "AliITSRecPoint.h"
14
15 Int_t AliITSStoreFindableTracksCompiled
16 (Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0)
17 {
18         // Make sure that ALICE objects are loaded
19         if (gAlice) {
20                 delete gAlice;
21                 gAlice = 0;
22         }
23         
24         // Define the names of all involved files
25         TString strEventFile(evname);
26         TString strOutputFile(evname);
27         strEventFile.Append(".root");
28         strOutputFile.Append("_tracks_");
29         strOutputFile += evnum;
30         strOutputFile.Append(".root");
31                 
32         // Connect the Root Galice file containing Geometry, Kine and Hits
33         TFile *fileEvent = (TFile*)gROOT->GetListOfFiles()->FindObject(strEventFile);
34         if (!fileEvent) fileEvent = new TFile(strEventFile,"UPDATE");
35         
36         // Get AliRun object from file
37         gAlice = (AliRun*)fileEvent->Get("gAlice");
38         if (gAlice) cout << "OK, found an AliRun object in file" << endl;
39         
40         // Get ITS related objects and data
41         AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS");
42         if (!ITS) {
43                 cerr << "ITS object not found!" << endl;
44                 return 1;
45         }
46         AliITSgeom *geometry = ITS->GetITSgeom(); 
47         if (!geometry) {
48                 cerr << "ITS geometry object not found!" << endl;
49                 return 2;
50         }
51         
52         // Count the number of modules per layer
53         Int_t nLadders, nDetectors, mod_min[6], mod_max[6];
54         for(Int_t i = 0; i < 6; i++) {
55                 nLadders = geometry->GetNladders(i + 1);
56                 nDetectors = geometry->GetNdetectors(i + 1);
57                 mod_min[i] = geometry->GetModuleIndex(i + 1, 1, 1);
58                 mod_max[i] = geometry->GetModuleIndex(i + 1, nLadders, nDetectors);
59         }
60         
61         // Load event and ITS recpoints
62         Int_t nParticles = gAlice->GetEvent(evnum);
63         cout << "Event number: " << evnum << endl;
64         cout << "# particles : " << nParticles <<endl;
65         if (nParticles <= 0) {
66                 cerr << "Can't have <= 0 particles!" << endl;
67                 return 3; 
68         }
69         AliITSRecPoint *recp = 0;
70         TClonesArray  *recPoints = ITS->RecPoints();
71         TObjArray     *particles = gAlice->Particles();
72         Int_t nTracks = gAlice->GetNtrack(); //FCA correction
73         Bool_t *hitITSLayer[6];
74         for (Int_t i = 0; i < 6; i++) {
75                 hitITSLayer[i] = new Bool_t[nTracks];
76                 for (Int_t j = 0; j < nTracks; j++) hitITSLayer[i][j] = kFALSE;
77         }
78         
79         // Load recpoints in event
80         TTree *TR = gAlice->TreeR();
81         if (!TR) {
82                 cerr << "TreeR object not found!" << endl;
83                 return 4;
84         }
85         
86         // Scan recpoints and define findable tracks
87         Int_t nModules = (Int_t)TR->GetEntries(), nPoints = 0;
88         cout << "Found " << nModules;
89         cout << " entries in the TreeR (must be one per module!)" << endl;
90         for (Int_t layer = 1; layer <= 6; layer++) {
91                 for (Int_t mod = mod_min[layer - 1]; mod <= mod_max[layer - 1]; mod++) {
92                         ITS->ResetRecPoints();
93                         TR->GetEntry(mod);
94                         nPoints = recPoints->GetEntries();
95                         if(!nPoints) {
96                                 cout << "Module " << mod << " is empty..." << endl;
97                                 continue;
98                         }
99                         for (Int_t point = 0; point < nPoints; point++) {
100                                 recp = (AliITSRecPoint*)recPoints->UncheckedAt(point);
101                                 for (Int_t it = 0; it < 3; it++) {
102                                         Int_t track = recp->GetLabel(it);
103                                         if(track < 0) continue;
104                                         if(track > nTracks) {
105                                                 cout << "Found track index " << track;
106                                                 cout << " whilw gAlice->GetNtrack() = " << nTracks << endl;
107                                                 continue;
108                                         }
109                                         hitITSLayer[layer - 1][track] = kTRUE;
110                                 } // loop over recpoint labels
111                         } //loop over points
112                 } //loop over modules
113         } //loop over layers
114         
115         // Scan the file of tracks in TPC to retrieve the findable TPC tracks
116         TString strLabelsTPC;
117         Int_t label, pdg_code, nFindablesTPC = 0;
118         Double_t dummy;
119         ifstream tpc("good_tracks_tpc");
120         while (tpc >> label >> pdg_code) {
121                 for (Int_t i = 0; i < 6; i++) tpc >> dummy;
122                 nFindablesTPC++;
123                 strLabelsTPC.Append(Form("[%d]", label));               
124         }
125                         
126         // Define the TTree with tracks data by means of a set of variables
127         Int_t nFindablesITS = 0, nFindablesITSTPC = 0;
128         Int_t nhits, tpc_ok, entry = 0;
129         Double_t vx, vy, vz;
130         Double_t px, py, pz, pt;
131
132         TTree *tree = new TTree("Tracks", "Findable tracks in ITS");
133         
134         tree->Branch("vx", &vx, "vx/D");
135         tree->Branch("vy", &vy, "vy/D");
136         tree->Branch("vz", &vz, "vz/D");
137         tree->Branch("px", &px, "px/D");
138         tree->Branch("py", &py, "py/D");
139         tree->Branch("pz", &pz, "pz/D");
140         tree->Branch("pt", &pt, "pt/D");
141         tree->Branch("label", &label, "label/I");
142         tree->Branch("entry", &entry, "entry/I");
143         tree->Branch("pdg_code", &pdg_code, "pdg_code/I");
144         tree->Branch("nhits", &nhits, "nhits/I");
145         tree->Branch("tpc_ok", &tpc_ok, "tpc_ok/I");
146         
147         // Fill the tree
148         cout << endl;
149         TParticle *p = 0;
150         for (Int_t i = 0; i < nTracks; i++) {
151                 p = gAlice->Particle(i);
152                 px = p->Px();
153                 py = p->Py();
154                 pz = p->Pz();
155                 pt = p->Pt();
156                 vx = p->Vx();
157                 vy = p->Vy();
158                 vz = p->Vz();
159                 nhits = 0;
160                 for (Int_t j = 0; j < 6; j++) if (hitITSLayer[j][i]) nhits++;
161                 if (nhits < nMinClusters) continue;
162                 cout << "Track " << i << " stored\r" << flush;
163                 tpc_ok = (strLabelsTPC.Contains(Form("[%d]", i)));
164                 pdg_code = p->GetPdgCode();
165                 label = i;
166                 nFindablesITS++;
167                 if (tpc_ok) nFindablesITSTPC++;
168                 tree->Fill();
169                 entry++;
170         }
171         
172         // Save into a file
173         TFile *fileOutput = new TFile(strOutputFile, "recreate");
174         tree->Write();
175         fileOutput->Close();
176         
177         cout << "# findable tracks in TPC     : " << nFindablesTPC << endl;
178         cout << "# findable tracks in ITS     : " << nFindablesITS << endl;
179         cout << "# findable tracks in ITS+TPC : " << nFindablesITSTPC << endl;
180         
181         return 0;
182 }