Update for Ds
[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"
d5eae813 8#include "AliITSRecPointContainer.h"
a3a3a28f 9#include "AliITSLoader.h"
10#include "AliLog.h"
11#include "AliRawReader.h"
12#include "AliRawReaderDate.h"
13#include "AliRawReaderRoot.h"
14#include "AliRunLoader.h"
15#include "AliITSVertexer3D.h"
37add1d1 16#include "AliITSVertexer3DTapan.h"
a3a3a28f 17#include "AliESDVertex.h"
18#include "AliMeanVertex.h"
a3a3a28f 19
199775b9 20const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
a3a3a28f 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//______________________________________________________________________
d5eae813 42AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
37add1d1 43fDetTypeRec(NULL),
44fVertexXY(NULL),
45fVertexZ(NULL),
a3a3a28f 46fNoEventsContr(0),
199775b9 47fTotContributors(0.),
48fAverContributors(0.),
a3a3a28f 49fFilterOnContributors(0),
d5eae813 50fMode(mode),
199775b9 51fVertexer(NULL),
52fAccEvents(fgkMaxNumOfEvents),
53fVertArray("AliESDVertex",fgkMaxNumOfEvents),
54fClu0(NULL),
55fIndex(0),
56fErrXCut(0.),
57fRCut(0.),
58fLowSPD0(0),
59fHighSPD0(0)
a3a3a28f 60{
61 // Default Constructor
199775b9 62 Reset();
a3a3a28f 63
37add1d1 64 // Histograms initialization
8d7f6205 65 const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
37add1d1 66 const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
67 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
68 2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
69 2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
63845149 70 fVertexXY->SetXTitle("X , cm");
71 fVertexXY->SetYTitle("Y , cm");
72 fVertexXY->SetOption("colz");
37add1d1 73 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
74 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
63845149 75 fVertexZ->SetXTitle("Z , cm");
199775b9 76 fClu0 = new UInt_t [fgkMaxNumOfEvents];
77 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
78 fAccEvents.ResetAllBits(kTRUE);
a3a3a28f 79}
80
81//______________________________________________________________________
199775b9 82void AliITSMeanVertexer::Reset(){
83 // Reset averages
84 for(Int_t i=0;i<3;i++){
85 fWeighPosSum[i] = 0.;
86 fWeighSigSum[i] = 0.;
87 fAverPosSum[i] = 0.;
88 fWeighPos[i] = 0.;
89 fWeighSig[i] = 0.;
90 fAverPos[i] = 0.;
91 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
92 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
93 fNoEventsContr=0;
94 fTotContributors = 0.;
95 if(fVertexXY)fVertexXY->Reset();
96 if(fVertexZ)fVertexZ->Reset();
97 }
98}
99//______________________________________________________________________
37add1d1 100Bool_t AliITSMeanVertexer::Init() {
101 // Initialize filters
a3a3a28f 102 // Initialize geometry
37add1d1 103 // Initialize ITS classes
a3a3a28f 104
37add1d1 105 AliGeomManager::LoadGeometry();
106 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f 107
108 AliITSInitGeometry initgeom;
109 AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1 110 if (!geom) return kFALSE;
a3a3a28f 111 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1 112
cf0b1cea 113 fDetTypeRec = new AliITSDetTypeRec();
762133f4 114 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
37add1d1 115 fDetTypeRec->SetITSgeom(geom);
116 fDetTypeRec->SetDefaults();
117 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
118
a3a3a28f 119 // Initialize filter values to their defaults
120 SetFilterOnContributors();
199775b9 121 SetCutOnErrX();
122 SetCutOnR();
123 SetCutOnCls();
a3a3a28f 124
d5eae813 125 // Instatiate vertexer
126 if (!fMode) {
127 fVertexer = new AliITSVertexer3DTapan(1000);
128 }
129 else {
130 fVertexer = new AliITSVertexer3D();
131 fVertexer->SetDetTypeRec(fDetTypeRec);
199775b9 132 AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
133 alias->SetWideFiducialRegion(40.,0.5);
134 alias->SetNarrowFiducialRegion(0.5,0.5);
135 alias->SetDeltaPhiCuts(0.5,0.025);
136 alias->SetDCACut(0.1);
137 alias->SetPileupAlgo(3);
138 fVertexer->SetComputeMultiplicity(kFALSE);
d5eae813 139 }
37add1d1 140 return kTRUE;
a3a3a28f 141}
142
143//______________________________________________________________________
144AliITSMeanVertexer::~AliITSMeanVertexer() {
145 // Destructor
37add1d1 146 delete fDetTypeRec;
147 delete fVertexXY;
148 delete fVertexZ;
d5eae813 149 delete fVertexer;
a3a3a28f 150}
151
152//______________________________________________________________________
d5eae813 153Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
a3a3a28f 154 // Performs SPD local reconstruction
37add1d1 155 // and vertex finding
156 // returns true in case a vertex is found
157
158 // Run SPD cluster finder
d5eae813 159 AliITSRecPointContainer::Instance()->PrepareToRead();
37add1d1 160 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
161 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
162
163 Bool_t vtxOK = kFALSE;
199775b9 164 UInt_t nl1=0;
d5eae813 165 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
166 if (!fMode) {
37add1d1 167 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
a3a3a28f 168 }
37add1d1 169 else {
199775b9 170 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
171 nl1=rpcont->GetNClustersInLayerFast(1);
172 if(Filter(vtx,nl1)) vtxOK = kTRUE;
a3a3a28f 173 }
37add1d1 174 delete clustersTree;
199775b9 175 if (vtxOK && fMode){
176 new(fVertArray[fIndex])AliESDVertex(*vtx);
177 fClu0[fIndex]=nl1;
178 fIndex++;
179 }
37add1d1 180 if (vtx) delete vtx;
a3a3a28f 181
37add1d1 182 return vtxOK;
a3a3a28f 183}
184
185//______________________________________________________________________
37add1d1 186void AliITSMeanVertexer::WriteVertices(const char *filename){
187 // Compute mean vertex and
188 // store it along with the histograms
189 // in a file
190
191 TFile fmv(filename,"update");
199775b9 192 Bool_t itisOK = kFALSE;
193 if(ComputeMean(kFALSE)){
194 if(ComputeMean(kTRUE)){
195 Double_t cov[6];
196 cov[0] = fAverPosSq[0][0]; // variance x
197 cov[1] = fAverPosSq[0][1]; // cov xy
198 cov[2] = fAverPosSq[1][1]; // variance y
199 cov[3] = fAverPosSq[0][2]; // cov xz
200 cov[4] = fAverPosSq[1][2]; // cov yz
201 cov[5] = fAverPosSq[2][2]; // variance z
202 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
203 mv.SetTitle("Mean Vertex");
204 mv.SetName("MeanVertex");
205 // we have to add chi2 here
206 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
37add1d1 207
199775b9 208 mv.Write(mv.GetName(),TObject::kOverwrite);
209 vtx.Write(vtx.GetName(),TObject::kOverwrite);
210 itisOK = kTRUE;
211 }
a3a3a28f 212 }
199775b9 213
214 if(!itisOK){
a3a3a28f 215 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
216 }
37add1d1 217
218 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
219 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
220 fmv.Close();
a3a3a28f 221}
222
223//______________________________________________________________________
199775b9 224Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
a3a3a28f 225 // Apply selection criteria to events
226 Bool_t status = kFALSE;
199775b9 227 if(!vert)return status;
8d7f6205 228 // Remove vertices reconstructed with vertexerZ
229 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
a3a3a28f 230 Int_t ncontr = vert->GetNContributors();
a3a3a28f 231 AliDebug(1,Form("Number of contributors = %d",ncontr));
199775b9 232 Double_t ex = vert->GetXRes();
233 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
234 ((ex*1000.)<fErrXCut)) {
235 status = kTRUE;
236 fAccEvents.SetBitNumber(fIndex);
237 }
a3a3a28f 238 return status;
239}
240
241//______________________________________________________________________
242void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
243 // update mean vertex
244 Double_t currentPos[3],currentSigma[3];
245 vert->GetXYZ(currentPos);
246 vert->GetSigmaXYZ(currentSigma);
247 Bool_t goon = kTRUE;
248 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
249 if(!goon)return;
250 for(Int_t i=0;i<3;i++){
762133f4 251 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
252 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
253 fAverPosSum[i]+=currentPos[i];
a3a3a28f 254 }
255 for(Int_t i=0;i<3;i++){
256 for(Int_t j=i;j<3;j++){
762133f4 257 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f 258 }
259 }
37add1d1 260
199775b9 261 fTotContributors+=vert->GetNContributors();
37add1d1 262 fVertexXY->Fill(currentPos[0],currentPos[1]);
263 fVertexZ->Fill(currentPos[2]);
264
a3a3a28f 265 fNoEventsContr++;
266}
267
268//______________________________________________________________________
199775b9 269 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
270 // compute mean vertex
271 // called once with killOutliers = kFALSE and then again with
272 // killOutliers = kTRUE
273 Double_t wpos[3];
274 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
275 Reset();
276 for(Int_t i=0;i<fVertArray.GetEntries();i++){
277 AliESDVertex* vert = NULL;
762133f4 278
199775b9 279 // second pass
280 if(killOutliers && fAccEvents.TestBitNumber(i)){
281 vert=(AliESDVertex*)fVertArray[i];
282 Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
283 dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
284 dist=sqrt(dist)*10.; // distance in mm
285 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
286 }
287
288 if(!fAccEvents.TestBitNumber(i))continue;
289 vert=(AliESDVertex*)fVertArray[i];
290 AddToMean(vert);
291 }
292 if(fNoEventsContr < 2) return kFALSE;
293 Double_t nevents = fNoEventsContr;
294 for(Int_t i=0;i<3;i++){
295 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
296 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
297 fAverPos[i] = fAverPosSum[i]/nevents;
298 }
299 for(Int_t i=0;i<3;i++){
300 for(Int_t j=i;j<3;j++){
301 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
302 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
303 }
304 }
305
306 fAverContributors = fTotContributors/nevents;
307 return kTRUE;
308 }