Removal of the run-loaders from the algorithmic part of the vertexers code. Some...
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
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();
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     TTree* cltree = loader->TreeR();
253     vert = dovert->FindVertexForCurrentEvent(cltree);
254     AliMultiplicity *mult = dovert->GetMultiplicity();
255     if(Filter(vert,mult)){
256       AddToMean(vert); 
257       if(fWriteVertices){
258         loader->PostVertex(vert);
259         loader->WriteVertices();
260       }
261     }
262     else {
263       if(vert)delete vert;
264     }
265   }
266   if(ComputeMean()){
267     Double_t cov[6];
268     cov[0] =  fAverPosSq[0][0];  // variance x
269     cov[1] =  fAverPosSq[0][1];  // cov xy
270     cov[2] =  fAverPosSq[1][1];  // variance y
271     cov[3] =  fAverPosSq[0][2];  // cov xz
272     cov[4] =  fAverPosSq[1][2];  // cov yz
273     cov[5] =  fAverPosSq[2][2];  // variance z
274     AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
275     mv.SetTitle("Mean Vertex");
276     mv.SetName("Meanvertex");
277     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
278     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
279     TFile fmv(fMVFileName.Data(),"recreate");
280     mv.Write();
281     fmv.Close();
282   }
283   else {
284     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
285   }
286   delete dovert;
287 }
288
289 //______________________________________________________________________
290 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
291   // Apply selection criteria to events
292   Bool_t status = kFALSE;
293   if(!vert || !mult)return status;
294   Int_t ncontr = vert->GetNContributors();
295   Int_t ntracklets = mult->GetNumberOfTracklets();
296   AliDebug(1,Form("Number of contributors = %d",ncontr));
297   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
298   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
299   fAverTracklets += ntracklets;
300   fSigmaOnAverTracks += ntracklets*ntracklets;
301   return status;
302 }
303
304 //______________________________________________________________________
305 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
306   // update mean vertex
307   Double_t currentPos[3],currentSigma[3];
308   vert->GetXYZ(currentPos);
309   vert->GetSigmaXYZ(currentSigma);
310   Bool_t goon = kTRUE;
311   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
312   if(!goon)return;
313   for(Int_t i=0;i<3;i++){
314     fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
315     fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
316     fAverPos[i]+=currentPos[i];
317   }
318   for(Int_t i=0;i<3;i++){
319     for(Int_t j=i;j<3;j++){
320       fAverPosSq[i][j] += currentPos[i] * currentPos[j];
321     }
322   }
323   fNoEventsContr++;
324 }
325
326 //______________________________________________________________________
327 Bool_t AliITSMeanVertexer::ComputeMean(){
328   // compute mean vertex
329   if(fNoEventsContr < 2) return kFALSE;
330   Double_t nevents = fNoEventsContr;
331   for(Int_t i=0;i<3;i++){
332     fWeighPos[i] /= fWeighSig[i]; 
333     fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
334     fAverPos[i] /= nevents;
335   } 
336   for(Int_t i=0;i<3;i++){
337     for(Int_t j=i;j<3;j++){
338       fAverPosSq[i][j] /= (nevents -1.);
339       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
340     }
341   } 
342   fTotTracklets = fAverTracklets;  //  total number of tracklets used 
343   fAverTracklets /= nevents;
344   fSigmaOnAverTracks /= (nevents - 1);
345   Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
346   fSigmaOnAverTracks -= tmp;
347   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
348   return kTRUE;
349 }