Diamond for vertexer Z set in vertexer3D to the value defined in AliITSRecoParam...
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
CommitLineData
a3a3a28f 1#include <TFile.h>
37add1d1 2#include <TH1.h>
3#include <TH2.h>
a3a3a28f 4#include "AliGeomManager.h"
a3a3a28f 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"
37add1d1 15#include "AliITSVertexer3DTapan.h"
a3a3a28f 16#include "AliESDVertex.h"
17#include "AliMeanVertex.h"
18#include "AliMultiplicity.h"
19
a3a3a28f 20ClassImp(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//______________________________________________________________________
41AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
37add1d1 42fDetTypeRec(NULL),
43fVertexXY(NULL),
44fVertexZ(NULL),
a3a3a28f 45fNoEventsContr(0),
46fTotTracklets(0.),
47fAverTracklets(0.),
48fSigmaOnAverTracks(0.),
49fFilterOnContributors(0),
37add1d1 50fFilterOnTracklets(0)
a3a3a28f 51{
52 // Default Constructor
53 for(Int_t i=0;i<3;i++){
54 fWeighPos[i] = 0.;
55 fWeighSig[i] = 0.;
56 fAverPos[i] = 0.;
57 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
58 }
a3a3a28f 59
37add1d1 60 // Histograms initialization
61 const Float_t xLimit = 2.0, yLimit = 2.0, zLimit = 50.0;
62 const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
63 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
64 2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
65 2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
66 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
67 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
a3a3a28f 68}
69
70//______________________________________________________________________
37add1d1 71Bool_t AliITSMeanVertexer::Init() {
72 // Initialize filters
a3a3a28f 73 // Initialize geometry
37add1d1 74 // Initialize ITS classes
a3a3a28f 75
37add1d1 76 AliGeomManager::LoadGeometry();
77 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f 78
79 AliITSInitGeometry initgeom;
80 AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1 81 if (!geom) return kFALSE;
a3a3a28f 82 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1 83
cf0b1cea 84 fDetTypeRec = new AliITSDetTypeRec();
37add1d1 85 fDetTypeRec->SetITSgeom(geom);
86 fDetTypeRec->SetDefaults();
87 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
88
a3a3a28f 89 // Initialize filter values to their defaults
90 SetFilterOnContributors();
91 SetFilterOnTracklets();
a3a3a28f 92
37add1d1 93 return kTRUE;
a3a3a28f 94}
95
96//______________________________________________________________________
97AliITSMeanVertexer::~AliITSMeanVertexer() {
98 // Destructor
37add1d1 99 delete fDetTypeRec;
100 delete fVertexXY;
101 delete fVertexZ;
a3a3a28f 102}
103
104//______________________________________________________________________
37add1d1 105Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
a3a3a28f 106 // Performs SPD local reconstruction
37add1d1 107 // and vertex finding
108 // returns true in case a vertex is found
109
110 // Run SPD cluster finder
111 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
112 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
113
114 Bool_t vtxOK = kFALSE;
115 AliESDVertex *vtx = NULL;
116 // Run Tapan's vertex-finder
117 if (!mode) {
118 AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
119 vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
120 delete vertexer1;
121 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
a3a3a28f 122 }
37add1d1 123 else {
124 // Run standard vertexer3d
125 AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
cf0b1cea 126 vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
37add1d1 127 AliMultiplicity *mult = vertexer2->GetMultiplicity();
128 delete vertexer2;
129 if(Filter(vtx,mult)) vtxOK = kTRUE;
a3a3a28f 130 }
37add1d1 131 delete clustersTree;
132 if (vtxOK) AddToMean(vtx);
133 if (vtx) delete vtx;
a3a3a28f 134
37add1d1 135 return vtxOK;
a3a3a28f 136}
137
138//______________________________________________________________________
37add1d1 139void AliITSMeanVertexer::WriteVertices(const char *filename){
140 // Compute mean vertex and
141 // store it along with the histograms
142 // in a file
143
144 TFile fmv(filename,"update");
145
a3a3a28f 146 if(ComputeMean()){
147 Double_t cov[6];
148 cov[0] = fAverPosSq[0][0]; // variance x
149 cov[1] = fAverPosSq[0][1]; // cov xy
150 cov[2] = fAverPosSq[1][1]; // variance y
151 cov[3] = fAverPosSq[0][2]; // cov xz
152 cov[4] = fAverPosSq[1][2]; // cov yz
153 cov[5] = fAverPosSq[2][2]; // variance z
37add1d1 154 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
a3a3a28f 155 mv.SetTitle("Mean Vertex");
37add1d1 156 mv.SetName("MeanVertex");
a3a3a28f 157 AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
158 AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
37add1d1 159 // we have to add chi2 here
160 AliESDVertex vtx(fWeighPos,cov,0,fAverTracklets,"MeanVertexPos");
161
162 mv.Write(mv.GetName(),TObject::kOverwrite);
163 vtx.Write(vtx.GetName(),TObject::kOverwrite);
a3a3a28f 164 }
165 else {
166 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
167 }
37add1d1 168
169 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
170 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
171 fmv.Close();
a3a3a28f 172}
173
174//______________________________________________________________________
175Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
176 // Apply selection criteria to events
177 Bool_t status = kFALSE;
178 if(!vert || !mult)return status;
179 Int_t ncontr = vert->GetNContributors();
180 Int_t ntracklets = mult->GetNumberOfTracklets();
181 AliDebug(1,Form("Number of contributors = %d",ncontr));
182 AliDebug(1,Form("Number of tracklets = %d",ntracklets));
183 if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
184 fAverTracklets += ntracklets;
185 fSigmaOnAverTracks += ntracklets*ntracklets;
186 return status;
187}
188
189//______________________________________________________________________
190void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
191 // update mean vertex
192 Double_t currentPos[3],currentSigma[3];
193 vert->GetXYZ(currentPos);
194 vert->GetSigmaXYZ(currentSigma);
195 Bool_t goon = kTRUE;
196 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
197 if(!goon)return;
198 for(Int_t i=0;i<3;i++){
199 fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
200 fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
201 fAverPos[i]+=currentPos[i];
202 }
203 for(Int_t i=0;i<3;i++){
204 for(Int_t j=i;j<3;j++){
205 fAverPosSq[i][j] += currentPos[i] * currentPos[j];
206 }
207 }
37add1d1 208
209 fVertexXY->Fill(currentPos[0],currentPos[1]);
210 fVertexZ->Fill(currentPos[2]);
211
a3a3a28f 212 fNoEventsContr++;
213}
214
215//______________________________________________________________________
216Bool_t AliITSMeanVertexer::ComputeMean(){
217 // compute mean vertex
218 if(fNoEventsContr < 2) return kFALSE;
219 Double_t nevents = fNoEventsContr;
220 for(Int_t i=0;i<3;i++){
221 fWeighPos[i] /= fWeighSig[i];
222 fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
223 fAverPos[i] /= nevents;
224 }
225 for(Int_t i=0;i<3;i++){
226 for(Int_t j=i;j<3;j++){
227 fAverPosSq[i][j] /= (nevents -1.);
228 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
229 }
230 }
231 fTotTracklets = fAverTracklets; // total number of tracklets used
232 fAverTracklets /= nevents;
233 fSigmaOnAverTracks /= (nevents - 1);
234 Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
235 fSigmaOnAverTracks -= tmp;
236 fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
237 return kTRUE;
238}