]>
Commit | Line | Data |
---|---|---|
a3a3a28f | 1 | #include <TFile.h> |
2 | #include "AliGeomManager.h" | |
3 | #include "AliHeader.h" | |
4 | #include "AliITSDetTypeRec.h" | |
5 | #include "AliITSInitGeometry.h" | |
6 | #include "AliITSMeanVertexer.h" | |
7 | #include "AliITSLoader.h" | |
8 | #include "AliLog.h" | |
9 | #include "AliRawReader.h" | |
10 | #include "AliRawReaderDate.h" | |
11 | #include "AliRawReaderRoot.h" | |
12 | #include "AliRunLoader.h" | |
13 | #include "AliITSVertexer3D.h" | |
14 | #include "AliESDVertex.h" | |
15 | #include "AliMeanVertex.h" | |
16 | #include "AliMultiplicity.h" | |
17 | ||
18 | const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root"; | |
19 | ||
20 | ||
21 | ClassImp(AliITSMeanVertexer) | |
22 | ||
23 | /////////////////////////////////////////////////////////////////////// | |
24 | // // | |
25 | // Class to compute vertex position using SPD local reconstruction // | |
26 | // An average vertex position using all the events // | |
27 | // is built and saved // | |
28 | // Individual vertices can be optionally stored // | |
29 | // Origin: M.Masera (masera@to.infn.it) // | |
30 | // Usage: | |
31 | // AliITSMeanVertexer mv("RawDataFileName"); | |
32 | // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root | |
33 | // .... optional setters .... | |
34 | // mv.Reconstruct(); // SPD local reconstruction | |
35 | // mv.DoVertices(); | |
36 | // Resulting AliMeanVertex object is stored in a file named fMVFileName | |
37 | /////////////////////////////////////////////////////////////////////// | |
38 | ||
39 | /* $Id$ */ | |
40 | ||
41 | //______________________________________________________________________ | |
42 | AliITSMeanVertexer::AliITSMeanVertexer():TObject(), | |
43 | fLoaderFileName(), | |
44 | fGeometryFileName(), | |
45 | fMVFileName(), | |
46 | fRawReader(), | |
47 | fRunLoader(), | |
48 | fNoEventsContr(0), | |
49 | fTotTracklets(0.), | |
50 | fAverTracklets(0.), | |
51 | fSigmaOnAverTracks(0.), | |
52 | fFilterOnContributors(0), | |
53 | fFilterOnTracklets(0), | |
54 | fWriteVertices(kFALSE) | |
55 | { | |
56 | // Default Constructor | |
57 | for(Int_t i=0;i<3;i++){ | |
58 | fWeighPos[i] = 0.; | |
59 | fWeighSig[i] = 0.; | |
60 | fAverPos[i] = 0.; | |
61 | for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.; | |
62 | } | |
63 | ||
64 | } | |
65 | ||
66 | //______________________________________________________________________ | |
67 | AliITSMeanVertexer::AliITSMeanVertexer(TString filename):TObject(), | |
68 | fLoaderFileName(), | |
69 | fGeometryFileName(), | |
70 | fMVFileName(fgkMVFileNameDefault), | |
71 | fRawReader(), | |
72 | fRunLoader(), | |
73 | fNoEventsContr(0), | |
74 | fTotTracklets(0.), | |
75 | fAverTracklets(0.), | |
76 | fSigmaOnAverTracks(0.), | |
77 | fFilterOnContributors(0), | |
78 | fFilterOnTracklets(0), | |
79 | fWriteVertices(kTRUE) | |
80 | { | |
81 | // Standard constructor | |
82 | ||
83 | for(Int_t i=0;i<3;i++){ | |
84 | fWeighPos[i] = 0.; | |
85 | fWeighSig[i] = 0.; | |
86 | fAverPos[i] = 0.; | |
87 | for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.; | |
88 | } | |
89 | SetLoaderFileName(); | |
90 | SetGeometryFileName(); | |
91 | ||
92 | Init(filename); | |
93 | } | |
94 | //______________________________________________________________________ | |
95 | AliITSMeanVertexer::AliITSMeanVertexer(TString filename, | |
96 | TString loaderfilename, | |
97 | TString geometryfilename):TObject(), | |
98 | fLoaderFileName(), | |
99 | fGeometryFileName(), | |
100 | fMVFileName(fgkMVFileNameDefault), | |
101 | fRawReader(), | |
102 | fRunLoader(), | |
103 | fNoEventsContr(0), | |
104 | fTotTracklets(0.), | |
105 | fAverTracklets(0.), | |
106 | fSigmaOnAverTracks(0.), | |
107 | fFilterOnContributors(0), | |
108 | fFilterOnTracklets(0), | |
109 | fWriteVertices(kTRUE) | |
110 | { | |
111 | // Standard constructor with explicit geometry file name assignment | |
112 | for(Int_t i=0;i<3;i++){ | |
113 | fWeighPos[i] = 0.; | |
114 | fWeighSig[i] = 0.; | |
115 | fAverPos[i] = 0.; | |
116 | for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.; | |
117 | } | |
118 | SetLoaderFileName(loaderfilename); | |
119 | SetGeometryFileName(geometryfilename); | |
120 | Init(filename); | |
121 | } | |
122 | ||
123 | //______________________________________________________________________ | |
124 | void AliITSMeanVertexer::Init(TString filename){ | |
125 | // Initialization part common to different constructors | |
126 | if(filename.IsNull()){ | |
127 | AliFatal("Please, provide a valid file name for raw data file\n"); | |
128 | } | |
129 | // if file name ends with root a raw reader ROOT is assumed | |
130 | if(filename.EndsWith(".root")){ | |
131 | fRawReader = new AliRawReaderRoot(filename); | |
132 | } | |
133 | else { // DATE raw reader is assumed | |
134 | fRawReader = new AliRawReaderDate(filename); | |
135 | fRawReader->SelectEvents(7); | |
136 | } | |
137 | fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate"); | |
138 | fRunLoader->MakeTree("E"); | |
139 | Int_t iEvent = 0; | |
140 | while (fRawReader->NextEvent()) { | |
141 | fRunLoader->SetEventNumber(iEvent); | |
142 | fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), | |
143 | iEvent, iEvent); | |
144 | fRunLoader->MakeTree("H"); | |
145 | fRunLoader->TreeE()->Fill(); | |
146 | iEvent++; | |
147 | } | |
148 | fRawReader->RewindEvents(); | |
149 | fRunLoader->SetNumberOfEventsPerFile(iEvent); | |
150 | fRunLoader->WriteHeader("OVERWRITE"); | |
151 | Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS"); | |
152 | if(retval != 0)AliFatal("Not able to add ITS detector"); | |
153 | AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName()); | |
154 | fRunLoader->AddLoader(loader); | |
155 | fRunLoader->CdGAFile(); | |
156 | fRunLoader->Write(0, TObject::kOverwrite); | |
157 | // Initialize geometry | |
158 | ||
159 | AliGeomManager::LoadGeometry(fGeometryFileName.Data()); | |
160 | ||
161 | AliITSInitGeometry initgeom; | |
162 | AliITSgeom *geom = initgeom.CreateAliITSgeom(); | |
163 | printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data()); | |
164 | loader->SetITSgeom(geom); | |
165 | // Initialize filter values to their defaults | |
166 | SetFilterOnContributors(); | |
167 | SetFilterOnTracklets(); | |
168 | } | |
169 | ||
170 | //______________________________________________________________________ | |
171 | AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr), | |
172 | fLoaderFileName(vtxr.fLoaderFileName), | |
173 | fGeometryFileName(vtxr.fGeometryFileName), | |
174 | fMVFileName(vtxr.fMVFileName), | |
175 | fRawReader(vtxr.fRawReader), | |
176 | fRunLoader(vtxr.fRunLoader), | |
177 | fNoEventsContr(vtxr.fNoEventsContr), | |
178 | fTotTracklets(vtxr.fTotTracklets), | |
179 | fAverTracklets(vtxr.fAverTracklets), | |
180 | fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks), | |
181 | fFilterOnContributors(vtxr.fFilterOnContributors), | |
182 | fFilterOnTracklets(vtxr.fFilterOnTracklets), | |
183 | fWriteVertices(vtxr.fWriteVertices) | |
184 | { | |
185 | // Copy constructor | |
186 | // Copies are not allowed. The method is protected to avoid misuse. | |
187 | AliFatal("Copy constructor not allowed\n"); | |
188 | ||
189 | } | |
190 | ||
191 | //______________________________________________________________________ | |
192 | AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){ | |
193 | // Assignment operator | |
194 | AliError("Assignment operator not allowed\n"); | |
195 | return *this; | |
196 | } | |
197 | ||
198 | //______________________________________________________________________ | |
199 | AliITSMeanVertexer::~AliITSMeanVertexer() { | |
200 | // Destructor | |
201 | delete fRawReader; | |
202 | delete fRunLoader; | |
203 | ||
204 | } | |
205 | ||
206 | //______________________________________________________________________ | |
207 | void AliITSMeanVertexer::Reconstruct(){ | |
208 | // Performs SPD local reconstruction | |
209 | AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader")); | |
210 | if (!loader) { | |
211 | AliFatal("ITS loader not found"); | |
212 | return; | |
213 | } | |
214 | AliITSDetTypeRec* rec = new AliITSDetTypeRec(); | |
215 | rec->SetITSgeom(loader->GetITSgeom()); | |
216 | rec->SetDefaults(); | |
217 | ||
218 | rec->SetDefaultClusterFindersV2(kTRUE); | |
219 | ||
220 | Int_t nEvents = fRunLoader->GetNumberOfEvents(); | |
221 | fRawReader->RewindEvents(); | |
222 | for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) { | |
223 | fRawReader->NextEvent(); | |
224 | fRunLoader->GetEvent(iEvent); | |
225 | AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent)); | |
226 | loader->LoadRecPoints("update"); | |
227 | loader->CleanRecPoints(); | |
228 | loader->MakeRecPointsContainer(); | |
229 | TTree *tR = loader->TreeR(); | |
230 | if(!tR){ | |
231 | AliFatal("Tree R pointer not found - Abort \n"); | |
232 | break; | |
233 | } | |
234 | rec->DigitsToRecPoints(fRawReader,tR,"SPD"); | |
235 | rec->ResetRecPoints(); | |
236 | rec->ResetClusters(); | |
237 | loader->WriteRecPoints("OVERWRITE"); | |
238 | loader->UnloadRecPoints(); | |
239 | } | |
240 | ||
241 | } | |
242 | ||
243 | //______________________________________________________________________ | |
244 | void AliITSMeanVertexer::DoVertices(){ | |
245 | // Loop on all events and compute 3D vertices | |
246 | AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader")); | |
247 | AliITSVertexer3D *dovert = new AliITSVertexer3D("ITS.Vert3D.root"); | |
248 | AliESDVertex *vert = 0; | |
249 | Int_t nevents = fRunLoader->TreeE()->GetEntries(); | |
250 | for(Int_t i=0; i<nevents; i++){ | |
251 | fRunLoader->GetEvent(i); | |
252 | vert = dovert->FindVertexForCurrentEvent(i); | |
253 | AliMultiplicity *mult = dovert->GetMultiplicity(); | |
254 | if(Filter(vert,mult)){ | |
255 | AddToMean(vert); | |
256 | if(fWriteVertices){ | |
257 | loader->PostVertex(vert); | |
258 | loader->WriteVertices(); | |
259 | } | |
260 | } | |
261 | else { | |
262 | if(vert)delete vert; | |
263 | } | |
264 | } | |
265 | if(ComputeMean()){ | |
266 | Double_t cov[6]; | |
267 | cov[0] = fAverPosSq[0][0]; // variance x | |
268 | cov[1] = fAverPosSq[0][1]; // cov xy | |
269 | cov[2] = fAverPosSq[1][1]; // variance y | |
270 | cov[3] = fAverPosSq[0][2]; // cov xz | |
271 | cov[4] = fAverPosSq[1][2]; // cov yz | |
272 | cov[5] = fAverPosSq[2][2]; // variance z | |
273 | AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks); | |
274 | mv.SetTitle("Mean Vertex"); | |
275 | mv.SetName("Meanvertex"); | |
276 | AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets())); | |
277 | AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks())); | |
278 | TFile fmv(fMVFileName.Data(),"recreate"); | |
279 | mv.Write(); | |
280 | fmv.Close(); | |
281 | } | |
282 | else { | |
283 | AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr)); | |
284 | } | |
285 | delete dovert; | |
286 | } | |
287 | ||
288 | //______________________________________________________________________ | |
289 | Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){ | |
290 | // Apply selection criteria to events | |
291 | Bool_t status = kFALSE; | |
292 | if(!vert || !mult)return status; | |
293 | Int_t ncontr = vert->GetNContributors(); | |
294 | Int_t ntracklets = mult->GetNumberOfTracklets(); | |
295 | AliDebug(1,Form("Number of contributors = %d",ncontr)); | |
296 | AliDebug(1,Form("Number of tracklets = %d",ntracklets)); | |
297 | if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE; | |
298 | fAverTracklets += ntracklets; | |
299 | fSigmaOnAverTracks += ntracklets*ntracklets; | |
300 | return status; | |
301 | } | |
302 | ||
303 | //______________________________________________________________________ | |
304 | void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){ | |
305 | // update mean vertex | |
306 | Double_t currentPos[3],currentSigma[3]; | |
307 | vert->GetXYZ(currentPos); | |
308 | vert->GetSigmaXYZ(currentSigma); | |
309 | Bool_t goon = kTRUE; | |
310 | for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE; | |
311 | if(!goon)return; | |
312 | for(Int_t i=0;i<3;i++){ | |
313 | fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i]; | |
314 | fWeighSig[i]+=1./currentSigma[i]/currentSigma[i]; | |
315 | fAverPos[i]+=currentPos[i]; | |
316 | } | |
317 | for(Int_t i=0;i<3;i++){ | |
318 | for(Int_t j=i;j<3;j++){ | |
319 | fAverPosSq[i][j] += currentPos[i] * currentPos[j]; | |
320 | } | |
321 | } | |
322 | fNoEventsContr++; | |
323 | } | |
324 | ||
325 | //______________________________________________________________________ | |
326 | Bool_t AliITSMeanVertexer::ComputeMean(){ | |
327 | // compute mean vertex | |
328 | if(fNoEventsContr < 2) return kFALSE; | |
329 | Double_t nevents = fNoEventsContr; | |
330 | for(Int_t i=0;i<3;i++){ | |
331 | fWeighPos[i] /= fWeighSig[i]; | |
332 | fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]); | |
333 | fAverPos[i] /= nevents; | |
334 | } | |
335 | for(Int_t i=0;i<3;i++){ | |
336 | for(Int_t j=i;j<3;j++){ | |
337 | fAverPosSq[i][j] /= (nevents -1.); | |
338 | fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j]; | |
339 | } | |
340 | } | |
341 | fTotTracklets = fAverTracklets; // total number of tracklets used | |
342 | fAverTracklets /= nevents; | |
343 | fSigmaOnAverTracks /= (nevents - 1); | |
344 | Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets; | |
345 | fSigmaOnAverTracks -= tmp; | |
346 | fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks); | |
347 | return kTRUE; | |
348 | } |