Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
1 #include <TFile.h>
2 #include <TH1.h>
3 #include <TH2.h>
4 #include "AliGeomManager.h"
5 #include "AliITSDetTypeRec.h"
6 #include "AliITSInitGeometry.h"
7 #include "AliITSMeanVertexer.h"
8 #include "AliITSLoader.h"
9 #include "AliLog.h"
10 #include "AliRawReader.h"
11 #include "AliRawReaderDate.h"
12 #include "AliRawReaderRoot.h"
13 #include "AliRunLoader.h"
14 #include "AliITSVertexer3D.h"
15 #include "AliITSVertexer3DTapan.h"
16 #include "AliESDVertex.h"
17 #include "AliMeanVertex.h"
18 #include "AliMultiplicity.h"
19
20 ClassImp(AliITSMeanVertexer)
21
22 ///////////////////////////////////////////////////////////////////////
23 //                                                                   //
24 // Class to compute vertex position using SPD local reconstruction   //
25 // An average vertex position using all the events                   //
26 // is built and saved                                                //
27 // Individual vertices can be optionally stored                      //
28 // Origin: M.Masera  (masera@to.infn.it)                             //
29 // Usage:
30 // AliITSMeanVertexer mv("RawDataFileName");
31 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
32 // ....  optional setters ....
33 // mv.Reconstruct();  // SPD local reconstruction
34 // mv.DoVertices(); 
35 // Resulting AliMeanVertex object is stored in a file named fMVFileName
36 ///////////////////////////////////////////////////////////////////////
37
38 /* $Id$ */
39
40 //______________________________________________________________________
41 AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
42 fDetTypeRec(NULL),
43 fVertexXY(NULL),
44 fVertexZ(NULL),
45 fNoEventsContr(0),
46 fTotTracklets(0.),
47 fAverTracklets(0.),
48 fTotTrackletsSq(0.),
49 fSigmaOnAverTracks(0.), 
50 fFilterOnContributors(0),
51 fFilterOnTracklets(0)
52 {
53   // Default Constructor
54   for(Int_t i=0;i<3;i++){
55     fWeighPosSum[i] = 0.;
56     fWeighSigSum[i] = 0.;
57     fAverPosSum[i] = 0.;
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     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
63   }
64
65   // Histograms initialization
66   const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
67   const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
68   fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
69                        2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
70                        2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
71   fVertexXY->SetXTitle("X , cm");
72   fVertexXY->SetYTitle("Y , cm");
73   fVertexXY->SetOption("colz");
74   fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
75                        2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
76   fVertexZ->SetXTitle("Z , cm");
77 }
78
79 //______________________________________________________________________
80 Bool_t AliITSMeanVertexer::Init() {
81   // Initialize filters
82   // Initialize geometry
83   // Initialize ITS classes
84  
85   AliGeomManager::LoadGeometry();
86   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
87
88   AliITSInitGeometry initgeom;
89   AliITSgeom *geom = initgeom.CreateAliITSgeom();
90   if (!geom) return kFALSE;
91   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
92
93   fDetTypeRec = new AliITSDetTypeRec();
94   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
95   fDetTypeRec->SetITSgeom(geom);
96   fDetTypeRec->SetDefaults();
97   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
98
99   // Initialize filter values to their defaults
100   SetFilterOnContributors();
101   SetFilterOnTracklets();
102
103   return kTRUE;
104 }
105
106 //______________________________________________________________________
107 AliITSMeanVertexer::~AliITSMeanVertexer() {
108   // Destructor
109   delete fDetTypeRec;
110   delete fVertexXY;
111   delete fVertexZ;
112 }
113
114 //______________________________________________________________________
115 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
116   // Performs SPD local reconstruction
117   // and vertex finding
118   // returns true in case a vertex is found
119
120   // Run SPD cluster finder
121   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
122   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
123
124   Bool_t vtxOK = kFALSE;
125   AliESDVertex *vtx = NULL;
126   // Run Tapan's vertex-finder
127   if (!mode) {
128     AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
129     vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
130     delete vertexer1;
131     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
132   }
133   else {
134   // Run standard vertexer3d
135     AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
136     vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
137     AliMultiplicity *mult = vertexer2->GetMultiplicity();
138     delete vertexer2;
139     if(Filter(vtx,mult)) vtxOK = kTRUE;
140   }
141   delete clustersTree;
142   if (vtxOK) AddToMean(vtx);
143   if (vtx) delete vtx;
144
145   return vtxOK;
146 }
147
148 //______________________________________________________________________
149 void AliITSMeanVertexer::WriteVertices(const char *filename){
150   // Compute mean vertex and
151   // store it along with the histograms
152   // in a file
153   
154   TFile fmv(filename,"update");
155
156   if(ComputeMean()){
157     Double_t cov[6];
158     cov[0] =  fAverPosSq[0][0];  // variance x
159     cov[1] =  fAverPosSq[0][1];  // cov xy
160     cov[2] =  fAverPosSq[1][1];  // variance y
161     cov[3] =  fAverPosSq[0][2];  // cov xz
162     cov[4] =  fAverPosSq[1][2];  // cov yz
163     cov[5] =  fAverPosSq[2][2];  // variance z
164     AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
165     mv.SetTitle("Mean Vertex");
166     mv.SetName("MeanVertex");
167     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
168     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
169     // we have to add chi2 here
170     AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
171
172     mv.Write(mv.GetName(),TObject::kOverwrite);
173     vtx.Write(vtx.GetName(),TObject::kOverwrite);
174   }
175   else {
176     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
177   }
178
179   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
180   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
181   fmv.Close();
182 }
183
184 //______________________________________________________________________
185 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
186   // Apply selection criteria to events
187   Bool_t status = kFALSE;
188   if(!vert || !mult)return status;
189   // Remove vertices reconstructed with vertexerZ
190   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
191   Int_t ncontr = vert->GetNContributors();
192   Int_t ntracklets = mult->GetNumberOfTracklets();
193   AliDebug(1,Form("Number of contributors = %d",ncontr));
194   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
195   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
196   fTotTracklets += ntracklets;
197   fTotTrackletsSq += ntracklets*ntracklets;
198   return status;
199 }
200
201 //______________________________________________________________________
202 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
203   // update mean vertex
204   Double_t currentPos[3],currentSigma[3];
205   vert->GetXYZ(currentPos);
206   vert->GetSigmaXYZ(currentSigma);
207   Bool_t goon = kTRUE;
208   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
209   if(!goon)return;
210   for(Int_t i=0;i<3;i++){
211     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
212     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
213     fAverPosSum[i]+=currentPos[i];
214   }
215   for(Int_t i=0;i<3;i++){
216     for(Int_t j=i;j<3;j++){
217       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
218     }
219   }
220
221   fVertexXY->Fill(currentPos[0],currentPos[1]);
222   fVertexZ->Fill(currentPos[2]);
223
224   fNoEventsContr++;
225 }
226
227 //______________________________________________________________________
228 Bool_t AliITSMeanVertexer::ComputeMean(){
229   // compute mean vertex
230   if(fNoEventsContr < 2) return kFALSE;
231   Double_t nevents = fNoEventsContr;
232   for(Int_t i=0;i<3;i++){
233     fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; 
234     fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
235     fAverPos[i] = fAverPosSum[i]/nevents;
236   } 
237   for(Int_t i=0;i<3;i++){
238     for(Int_t j=i;j<3;j++){
239       fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
240       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
241     }
242   } 
243
244   fAverTracklets = fTotTracklets/nevents;
245   fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
246   fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
247   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
248   return kTRUE;
249 }