Replacing Header with Id
[u/mrichter/AliRoot.git] / STEER / AliTrackMapper.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////
19 //
20 // AliTrackMapper.cxx
21 // description: 
22 //        create an AliTrackMap - a relation between track label and 
23 //        it's index in TreeH. Check all branches with hits.
24 //        Output is in the root file.
25 //
26 ////////////////////////////////////////////////////////////////////////
27
28 #include <Riostream.h>
29
30 #include "TTree.h"
31 #include "TFile.h"
32 #include "TStopwatch.h"
33
34 #include "AliDetector.h"
35 #include "AliTrackMapper.h"
36 #include "AliTrackMap.h"
37 #include "AliRun.h"
38 #include "AliHit.h"
39
40 ClassImp(AliTrackMapper)
41
42 //_______________________________________________________________________
43 AliTrackMapper::AliTrackMapper(): 
44   fDEBUG(0)
45 {
46   //
47   // Default ctor
48   //
49 }
50
51
52 //_______________________________________________________________________
53 void AliTrackMapper::CreateMap(Int_t nEvents, Int_t firstEventNr,
54                                const char* fnMap, const char* fnHits)
55 {
56   //
57   // method to create a track map for a given number of events starting 
58   // from event firstEventNr. This method just opens and closes files and
59   // loops over events, the creation of map is delegated to 
60   // CreateMap(Int_t eventNr, TFile* fileMap) method
61   //
62   TStopwatch timer;
63   timer.Start();
64   
65   TFile *fileMap=TFile::Open(fnMap,"new");
66   if (!fileMap->IsOpen()) {cerr<<"Can't open output file "<<fnMap<<"!\n"; return;}
67   
68   TFile *fileHits=TFile::Open(fnHits);
69   if (!fileHits->IsOpen()) {cerr<<"Can't open input file "<<fnHits<<"!\n"; return;}
70   if (!(gAlice=dynamic_cast<AliRun*>(fileHits->Get("gAlice")))) {
71     cerr<<"gAlice have not been found on galice.root !\n";
72     return;
73   }
74   
75   for (Int_t eventNr = firstEventNr; eventNr < firstEventNr+nEvents;
76        eventNr++) {
77     CreateMap(eventNr,fileMap);
78   } // end loop over events
79   
80   delete gAlice;
81   gAlice = 0;
82   fileHits->Close();
83   delete fileHits;
84   fileMap->Close();
85   delete fileMap;
86   timer.Stop();
87   if (fDEBUG > 0) timer.Print();
88 }
89
90 //_______________________________________________________________________
91 Int_t  AliTrackMapper::CreateMap(Int_t eventNr, TFile* fileMap) const
92 {
93   //
94   // create an AliTrackMap for a given event
95   // correct gAlice must be already present in memory
96   //
97   Int_t nGenPrimPlusSecParticles = gAlice->GetEvent(eventNr);
98   if (fDEBUG > 1) cout<<"nGenPrimPlusSecParticles = "<<nGenPrimPlusSecParticles<<endl;
99   if (nGenPrimPlusSecParticles < 1) {
100     cerr<<"No primary particles found in event "<<eventNr<<endl;
101     return -1;
102   }
103
104   TTree *treeK = gAlice->TreeK();
105   if (!treeK) {
106     cerr<<"Error: Event "<<eventNr<<", treeK not found."<<endl;
107     return -1;
108   }
109   Int_t nAllParticles = static_cast<Int_t>(treeK->GetEntries());
110   Int_t *trackMap = new Int_t[nAllParticles];
111   for (Int_t i = 0; i<nAllParticles; i++) {trackMap[i] = kNoEntry;}
112
113   TTree *treeH = gAlice->TreeH();
114   if (!treeH) {
115     cerr<<"Error: Event "<<eventNr<<", treeH not found."<<endl;
116     return -1;
117   }
118   Int_t treeHEntries = static_cast<Int_t>(treeH->GetEntries());
119   if (fDEBUG > 1) cout<<"treeHEntries "<<treeHEntries<<endl;
120
121   TObjArray *modules = gAlice->Detectors();
122   if (!modules) {
123     cerr<<"TObjArray with modules not found."<<endl;
124     return -1;
125   }
126   Int_t nModules = static_cast<Int_t>(modules->GetEntries());
127   AliHit* hit;
128   for (Int_t treeHIndex = 0; treeHIndex < treeHEntries; treeHIndex++) {
129     gAlice->ResetHits();
130     treeH->GetEvent(treeHIndex);
131     for (Int_t iModule = 0; iModule < nModules; iModule++) {
132       AliDetector * detector = dynamic_cast<AliDetector*>(modules->At(iModule));
133       if (!detector) continue;
134       // process only detectors with shunt = 0
135       if (detector->GetIshunt()) continue; 
136
137       hit=dynamic_cast<AliHit*>(detector->FirstHit(-1));
138       Int_t lastLabel=-1, label;
139       for( ; hit; hit=dynamic_cast<AliHit*>(detector->NextHit()) ) {
140         label=hit->Track();     
141         if (lastLabel != label) {
142           if (label < 0 || label >= nAllParticles) {
143             cerr<<"Error: label out of range. ";
144             cerr<<"Event "<<eventNr<<" treeHIndex "<<treeHIndex<<" label = "<<label<<endl;
145             return -2;
146           }
147           if (trackMap[label] >=0 && trackMap[label] != treeHIndex) {
148             cerr<<"Error: different treeHIndex for label "<<label
149                 <<" indeces: "<<trackMap[label]<<" != "<<treeHIndex;
150             cerr<<" event "<<eventNr<<" detector "<<detector->GetName()<<endl;
151             return -3;
152           }
153           trackMap[label] = treeHIndex;
154           if (fDEBUG > 2) cout<<"treeHIndex, label = "<<treeHIndex<<" "<<label<<endl;
155           lastLabel = label;
156         }
157       }
158     }
159   }
160
161   if (fDEBUG > 2) {
162     for (Int_t i = 0; i < nAllParticles; i++) {
163       cout<<eventNr<<"\t"<<i<<"\t"<<trackMap[i]<<endl;
164     }
165   }
166   fileMap->cd();
167   AliTrackMap* trackMapObject = new AliTrackMap(nAllParticles, trackMap);
168   trackMapObject->SetEventNr(eventNr);
169   trackMapObject->Write();
170   delete trackMapObject;
171
172   delete [] trackMap;
173   return 0;
174 }
175   
176 //_______________________________________________________________________
177 AliTrackMap* AliTrackMapper::LoadTrackMap(Int_t eventNr, const char* fnMap, TFile* &fileMap) {
178   //
179   // read an AliTrackMap object for the given event eventNr from 
180   // the file fileMap
181   //
182   fileMap=TFile::Open(fnMap);
183   if (!fileMap->IsOpen()) {cerr<<"Can't open file "<<fnMap<<" with map!\n"; return 0;}
184   char mapName[20];
185   sprintf(mapName,"AliTrackMap_%5.5d",eventNr);
186   AliTrackMap* trackMapObject = dynamic_cast<AliTrackMap*>(fileMap->Get(mapName));
187   if (!trackMapObject) {
188     cerr<<"Error: map named "<<mapName<<" not found."<<endl;
189     return 0;
190   }
191   return trackMapObject;
192 }
193
194 //_______________________________________________________________________
195 void AliTrackMapper::CheckTrackMap(Int_t eventNr, const char* fnMap) {
196   //
197   // 
198   //
199   TFile *fileMap;
200   AliTrackMap* trackMapObject = LoadTrackMap(eventNr, fnMap, fileMap);
201   if (!trackMapObject) return;
202   
203   trackMapObject->PrintValues();
204   
205   delete trackMapObject;
206   fileMap->Close();
207   delete fileMap;
208 }
209