Processing SPD Mean Vertex only in PHYSICS runs.
[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 "TError.h"
31 #include "TFile.h"
32 #include "TStopwatch.h"
33 #include "TTree.h"
34
35 #include "AliLog.h"
36 #include "AliDetector.h"
37 #include "AliHit.h"
38 #include "AliLoader.h"
39 #include "AliRun.h"
40 #include "AliRunLoader.h"
41 #include "AliTrackMap.h"
42 #include "AliTrackMapper.h"
43
44 ClassImp(AliTrackMapper)
45
46 //_______________________________________________________________________
47 AliTrackMapper::AliTrackMapper(): 
48   fDEBUG(0)
49 {
50   //
51   // Default ctor
52   //
53 }
54
55
56 //_______________________________________________________________________
57 void AliTrackMapper::CreateMap(Int_t nEvents, Int_t firstEventNr,
58                                const char* fnMap, const char* fnHits)
59 {
60   //
61   // method to create a track map for a given number of events starting 
62   // from event firstEventNr. This method just opens and closes files and
63   // loops over events, the creation of map is delegated to 
64   // CreateMap(Int_t eventNr, TFile* fileMap) method
65   //
66   TStopwatch timer;
67   timer.Start();
68   
69   TFile *fileMap=TFile::Open(fnMap,"new");
70   if (!fileMap->IsOpen()) {AliErrorClass(Form("Can't open output file %s!", fnMap)); return;}
71
72   AliRunLoader* rl = AliRunLoader::Open(fnHits);
73   if (rl == 0x0) {AliErrorClass(Form("Can't open input file %s!\n", fnHits)); return;}
74   if (rl->LoadgAlice())
75     {
76       AliErrorClass("Error occured while loading AliRun");
77       return;
78     }
79     
80   if (!(gAlice=rl->GetAliRun())) {
81     AliErrorClass("gAlice have not been found on session !");
82     return;
83   }
84   
85   rl->LoadKinematics();
86   
87   for (Int_t eventNr = firstEventNr; eventNr < firstEventNr+nEvents;
88        eventNr++) {
89     CreateMap(eventNr,fileMap,rl);
90   } // end loop over events
91   
92   delete rl;
93   fileMap->Close();
94   delete fileMap;
95   timer.Stop();
96   if (fDEBUG > 0) timer.Print();
97 }
98
99 //_______________________________________________________________________
100 Int_t  AliTrackMapper::CreateMap(Int_t eventNr, TFile* fileMap,AliRunLoader* rl) 
101 {
102   //
103   // create an AliTrackMap for a given event
104   // correct gAlice must be already present in memory
105   //
106
107   rl->GetEvent(eventNr);
108
109   TTree *treeK = rl->TreeK();
110   if (!treeK) {
111     AliErrorClass(Form("Event %d, treeK not found.", eventNr));
112     return -1;
113   }
114   Int_t nAllParticles = static_cast<Int_t>(treeK->GetEntries());
115   Int_t *trackMap = new Int_t[nAllParticles];
116   for (Int_t i = 0; i<nAllParticles; i++) {trackMap[i] = kNoEntry;}
117
118
119   TObjArray *modules = gAlice->Detectors();
120   if (!modules) {
121     AliErrorClass("TObjArray with modules not found.");
122     return -1;
123   }
124   Int_t nModules = static_cast<Int_t>(modules->GetEntries());
125   AliHit* hit;
126   for (Int_t iModule = 0; iModule < nModules; iModule++) 
127    {
128     AliDetector * detector = dynamic_cast<AliDetector*>(modules->At(iModule));
129     if (!detector) continue;
130     AliLoader* loader = detector->GetLoader();
131     if (loader == 0x0)
132      {
133        AliWarningClass(Form("Can not get loader from detector %s.",
134                        detector->GetName()));
135        continue;
136      }
137     Int_t retval = loader->LoadHits();
138     if (retval) {
139       AliErrorClass(Form("Event %d: error occured while loading hits for %s",
140                          eventNr,detector->GetName()));
141       return -1;
142      }
143     
144     TTree *treeH = loader->TreeH();
145     if (!treeH) {
146       AliErrorClass(Form("Event %d: Can not get TreeH for %s",
147                          eventNr,detector->GetName()));
148       return -1;
149      }
150     Int_t treeHEntries = static_cast<Int_t>(treeH->GetEntries());
151     AliDebugClass(2, Form("treeHEntries %d", treeHEntries));
152      
153     detector->ResetHits();
154     
155     for (Int_t treeHIndex = 0; treeHIndex < treeHEntries; treeHIndex++)  
156      { // process only detectors with shunt = 0
157       treeH->GetEvent(treeHIndex);
158       if (detector->GetIshunt()) continue; 
159
160       hit=dynamic_cast<AliHit*>(detector->FirstHit(-1));
161       Int_t lastLabel=-1, label;
162       for( ; hit; hit=dynamic_cast<AliHit*>(detector->NextHit()) ) {
163        label=hit->Track();       
164        if (lastLabel != label) {
165          if (label < 0 || label >= nAllParticles) {
166            AliErrorClass(Form("label out of range. Event %d treeHIndex %d label = %d",
167                               eventNr, treeHIndex, label));
168            return -2;
169          }
170          if (trackMap[label] >=0 && trackMap[label] != treeHIndex) {
171            AliErrorClass(Form("different treeHIndex for label %d indeces: %d != %d event %d detector %s",
172                               label, trackMap[label], treeHIndex, eventNr, detector->GetName()));
173            return -3;
174          }
175          trackMap[label] = treeHIndex;
176          AliDebugClass(3, Form("treeHIndex, label = %d %d", treeHIndex, label));
177          lastLabel = label;
178        }
179       }
180     }//loop over hits in module
181     loader->UnloadHits();
182   }//loop over modules
183
184   ToAliDebugClass(3,
185     for (Int_t i = 0; i < nAllParticles; i++) {
186       cout<<eventNr<<"\t"<<i<<"\t"<<trackMap[i]<<endl;
187     }
188   )
189   fileMap->cd();
190   AliTrackMap* trackMapObject = new AliTrackMap(nAllParticles, trackMap);
191   trackMapObject->SetEventNr(eventNr);
192   trackMapObject->Write();
193   delete trackMapObject;
194
195   delete [] trackMap;
196   return 0;
197 }
198   
199 //_______________________________________________________________________
200 AliTrackMap* AliTrackMapper::LoadTrackMap(Int_t eventNr, const char* fnMap, TFile* &fileMap) {
201   //
202   // read an AliTrackMap object for the given event eventNr from 
203   // the file fileMap
204   //
205   fileMap=TFile::Open(fnMap);
206   if (!fileMap->IsOpen()) {AliErrorClass(Form("Can't open file %s with map!", fnMap)); return 0;}
207   char mapName[20];
208   sprintf(mapName,"AliTrackMap_%5.5d",eventNr);
209   AliTrackMap* trackMapObject = dynamic_cast<AliTrackMap*>(fileMap->Get(mapName));
210   if (!trackMapObject) {
211     AliErrorClass(Form("map named %s not found.", mapName));
212     return 0;
213   }
214   return trackMapObject;
215 }
216
217 //_______________________________________________________________________
218 void AliTrackMapper::CheckTrackMap(Int_t eventNr, const char* fnMap) {
219   //
220   // 
221   //
222   TFile *fileMap;
223   AliTrackMap* trackMapObject = LoadTrackMap(eventNr, fnMap, fileMap);
224   if (!trackMapObject) return;
225   
226   trackMapObject->PrintValues();
227   
228   delete trackMapObject;
229   fileMap->Close();
230   delete fileMap;
231 }
232