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