Remove wrong selction: lost in 27435
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
CommitLineData
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
18const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root";
19
20
21ClassImp(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//______________________________________________________________________
42AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
43fLoaderFileName(),
44fGeometryFileName(),
45fMVFileName(),
46fRawReader(),
47fRunLoader(),
48fNoEventsContr(0),
49fTotTracklets(0.),
50fAverTracklets(0.),
51fSigmaOnAverTracks(0.),
52fFilterOnContributors(0),
53fFilterOnTracklets(0),
54fWriteVertices(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//______________________________________________________________________
308c2f7c 67AliITSMeanVertexer::AliITSMeanVertexer(TString &filename):TObject(),
a3a3a28f 68fLoaderFileName(),
69fGeometryFileName(),
70fMVFileName(fgkMVFileNameDefault),
71fRawReader(),
72fRunLoader(),
73fNoEventsContr(0),
74fTotTracklets(0.),
75fAverTracklets(0.),
76fSigmaOnAverTracks(0.),
77fFilterOnContributors(0),
78fFilterOnTracklets(0),
79fWriteVertices(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//______________________________________________________________________
308c2f7c 95AliITSMeanVertexer::AliITSMeanVertexer(TString &filename,
96 TString &loaderfilename,
97 TString &geometryfilename):TObject(),
a3a3a28f 98fLoaderFileName(),
99fGeometryFileName(),
100fMVFileName(fgkMVFileNameDefault),
101fRawReader(),
102fRunLoader(),
103fNoEventsContr(0),
104fTotTracklets(0.),
105fAverTracklets(0.),
106fSigmaOnAverTracks(0.),
107fFilterOnContributors(0),
108fFilterOnTracklets(0),
109fWriteVertices(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//______________________________________________________________________
308c2f7c 124void AliITSMeanVertexer::Init(TString &filename){
a3a3a28f 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);
a3a3a28f 135 }
136 fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
137 fRunLoader->MakeTree("E");
138 Int_t iEvent = 0;
139 while (fRawReader->NextEvent()) {
140 fRunLoader->SetEventNumber(iEvent);
141 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
142 iEvent, iEvent);
143 fRunLoader->MakeTree("H");
144 fRunLoader->TreeE()->Fill();
145 iEvent++;
146 }
147 fRawReader->RewindEvents();
148 fRunLoader->SetNumberOfEventsPerFile(iEvent);
149 fRunLoader->WriteHeader("OVERWRITE");
150 Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS");
151 if(retval != 0)AliFatal("Not able to add ITS detector");
152 AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName());
153 fRunLoader->AddLoader(loader);
154 fRunLoader->CdGAFile();
155 fRunLoader->Write(0, TObject::kOverwrite);
156 // Initialize geometry
157
158 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
159
160 AliITSInitGeometry initgeom;
161 AliITSgeom *geom = initgeom.CreateAliITSgeom();
162 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
163 loader->SetITSgeom(geom);
164 // Initialize filter values to their defaults
165 SetFilterOnContributors();
166 SetFilterOnTracklets();
167}
168
169//______________________________________________________________________
170AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr),
171fLoaderFileName(vtxr.fLoaderFileName),
172fGeometryFileName(vtxr.fGeometryFileName),
173fMVFileName(vtxr.fMVFileName),
174fRawReader(vtxr.fRawReader),
175fRunLoader(vtxr.fRunLoader),
176fNoEventsContr(vtxr.fNoEventsContr),
177fTotTracklets(vtxr.fTotTracklets),
178fAverTracklets(vtxr.fAverTracklets),
179fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks),
180fFilterOnContributors(vtxr.fFilterOnContributors),
181fFilterOnTracklets(vtxr.fFilterOnTracklets),
182fWriteVertices(vtxr.fWriteVertices)
183{
184 // Copy constructor
185 // Copies are not allowed. The method is protected to avoid misuse.
186 AliFatal("Copy constructor not allowed\n");
187
188}
189
190//______________________________________________________________________
191AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){
192 // Assignment operator
193 AliError("Assignment operator not allowed\n");
194 return *this;
195}
196
197//______________________________________________________________________
198AliITSMeanVertexer::~AliITSMeanVertexer() {
199 // Destructor
200 delete fRawReader;
201 delete fRunLoader;
202
203}
204
205//______________________________________________________________________
206void AliITSMeanVertexer::Reconstruct(){
207 // Performs SPD local reconstruction
208 AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
209 if (!loader) {
210 AliFatal("ITS loader not found");
211 return;
212 }
213 AliITSDetTypeRec* rec = new AliITSDetTypeRec();
214 rec->SetITSgeom(loader->GetITSgeom());
215 rec->SetDefaults();
216
217 rec->SetDefaultClusterFindersV2(kTRUE);
218
219 Int_t nEvents = fRunLoader->GetNumberOfEvents();
220 fRawReader->RewindEvents();
221 for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
222 fRawReader->NextEvent();
223 fRunLoader->GetEvent(iEvent);
224 AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent));
225 loader->LoadRecPoints("update");
226 loader->CleanRecPoints();
227 loader->MakeRecPointsContainer();
228 TTree *tR = loader->TreeR();
229 if(!tR){
230 AliFatal("Tree R pointer not found - Abort \n");
231 break;
232 }
233 rec->DigitsToRecPoints(fRawReader,tR,"SPD");
234 rec->ResetRecPoints();
235 rec->ResetClusters();
236 loader->WriteRecPoints("OVERWRITE");
237 loader->UnloadRecPoints();
238 }
239
240}
241
242//______________________________________________________________________
243void AliITSMeanVertexer::DoVertices(){
244 // Loop on all events and compute 3D vertices
245 AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
308c2f7c 246 AliITSVertexer3D *dovert = new AliITSVertexer3D();
a3a3a28f 247 AliESDVertex *vert = 0;
248 Int_t nevents = fRunLoader->TreeE()->GetEntries();
249 for(Int_t i=0; i<nevents; i++){
250 fRunLoader->GetEvent(i);
308c2f7c 251 TTree* cltree = loader->TreeR();
252 vert = dovert->FindVertexForCurrentEvent(cltree);
a3a3a28f 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//______________________________________________________________________
289Bool_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//______________________________________________________________________
304void 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//______________________________________________________________________
326Bool_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}