fix compilation
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
CommitLineData
a3a3a28f 1#include <TFile.h>
37add1d1 2#include <TH1.h>
3#include <TH2.h>
c93255fe 4#include <TTree.h>
a3a3a28f 5#include "AliGeomManager.h"
a3a3a28f 6#include "AliITSDetTypeRec.h"
7#include "AliITSInitGeometry.h"
8#include "AliITSMeanVertexer.h"
d5eae813 9#include "AliITSRecPointContainer.h"
a3a3a28f 10#include "AliITSLoader.h"
11#include "AliLog.h"
12#include "AliRawReader.h"
13#include "AliRawReaderDate.h"
14#include "AliRawReaderRoot.h"
15#include "AliRunLoader.h"
16#include "AliITSVertexer3D.h"
37add1d1 17#include "AliITSVertexer3DTapan.h"
a3a3a28f 18#include "AliESDVertex.h"
19#include "AliMeanVertex.h"
8e2abaa0 20//#include "AliCodeTimer.h"
a3a3a28f 21
199775b9 22const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
a3a3a28f 23ClassImp(AliITSMeanVertexer)
a3a3a28f 24///////////////////////////////////////////////////////////////////////
25// //
26// Class to compute vertex position using SPD local reconstruction //
27// An average vertex position using all the events //
28// is built and saved //
29// Individual vertices can be optionally stored //
30// Origin: M.Masera (masera@to.infn.it) //
8e2abaa0 31// Usage: //
32// Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
a3a3a28f 33///////////////////////////////////////////////////////////////////////
a3a3a28f 34/* $Id$ */
a3a3a28f 35//______________________________________________________________________
d5eae813 36AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
37add1d1 37fDetTypeRec(NULL),
38fVertexXY(NULL),
39fVertexZ(NULL),
a3a3a28f 40fNoEventsContr(0),
199775b9 41fTotContributors(0.),
42fAverContributors(0.),
a3a3a28f 43fFilterOnContributors(0),
d5eae813 44fMode(mode),
199775b9 45fVertexer(NULL),
46fAccEvents(fgkMaxNumOfEvents),
47fVertArray("AliESDVertex",fgkMaxNumOfEvents),
48fClu0(NULL),
49fIndex(0),
50fErrXCut(0.),
51fRCut(0.),
4f40f207 52fZCutDiamond(0.),
199775b9 53fLowSPD0(0),
8e2abaa0 54fHighSPD0(0),
55fMultH(NULL),
56fErrXH(NULL),
57fMultHa(NULL),
58fErrXHa(NULL),
a1b1ad51 59fDistH(NULL),
60fContrH(NULL),
61fContrHa(NULL)
a3a3a28f 62{
63 // Default Constructor
4f40f207 64 SetZFiducialRegion(); // sets fZCutDiamond to its default value
a1b1ad51 65 Reset(kFALSE,kFALSE);
a3a3a28f 66
37add1d1 67 // Histograms initialization
4ed24269 68 const Float_t xLimit = 0.6, yLimit = 0.6, zLimit = 50.0;
69 const Float_t xDelta = 0.003, yDelta = 0.003, zDelta = 0.2;
37add1d1 70 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
71 2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
72 2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
63845149 73 fVertexXY->SetXTitle("X , cm");
74 fVertexXY->SetYTitle("Y , cm");
75 fVertexXY->SetOption("colz");
37add1d1 76 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
77 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
63845149 78 fVertexZ->SetXTitle("Z , cm");
51a579a9 79 fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
80 fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
81 fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
8e2abaa0 82 fErrXHa->SetLineColor(kRed);
51a579a9 83 fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
8e2abaa0 84 fMultHa->SetLineColor(kRed);
85 fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
a1b1ad51 86 fContrH = new TH1F("contrib","Number of contributors - before cuts",300,0.5,300.5);
87 fContrHa = new TH1F("contriba","Number of contributors - after cuts",300,0.5,300.5);
88 fContrHa->SetLineColor(kRed);
89
199775b9 90 fClu0 = new UInt_t [fgkMaxNumOfEvents];
91 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
92 fAccEvents.ResetAllBits(kTRUE);
a3a3a28f 93}
94
95//______________________________________________________________________
a1b1ad51 96void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
199775b9 97 // Reset averages
a1b1ad51 98 // Both arguments are kFALSE when the method is called by the constructor
99 // redefine2D is TRUE when the method is called by ComputeMean at the second
100 // pass.
101 // redefine2D is FALSE and complete is TRUE when the method is called
102 // after a failure cof ComputeMean(kTRUE)
103
4ed24269 104 static Int_t counter = 0;
a1b1ad51 105 if(fVertexXY && redefine2D){
4ed24269 106 counter++;
107 Double_t rangeX=fVertexXY->GetXaxis()->GetXmax()-fVertexXY->GetXaxis()->GetXmin();
108 Double_t rangeY=fVertexXY->GetYaxis()->GetXmax()-fVertexXY->GetYaxis()->GetXmin();
109 Int_t nx=fVertexXY->GetNbinsX();
110 Int_t ny=fVertexXY->GetNbinsY();
111 delete fVertexXY;
112 Double_t xmi=fWeighPos[0]-rangeX/2.;
113 Double_t xma=fWeighPos[0]+rangeX/2.;
114 Double_t ymi=fWeighPos[1]-rangeY/2.;
115 Double_t yma=fWeighPos[1]+rangeY/2.;
116 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",nx,xmi,xma,ny,ymi,yma);
117 fVertexXY->SetXTitle("X , cm");
118 fVertexXY->SetYTitle("Y , cm");
119 fVertexXY->SetOption("colz");
120 fVertexXY->SetDirectory(0);
121 }
a1b1ad51 122 else if(fVertexXY && !redefine2D){
123 fVertexXY->Reset();
124 }
125 if(fVertexZ){
126 fVertexZ->Reset();
127 fDistH->Reset();
128 if(complete){
129 fErrXH->Reset();
130 fMultH->Reset();
131 fErrXHa->Reset();
132 fMultHa->Reset();
133 fContrH->Reset();
134 fContrHa->Reset();
135 }
136 }
199775b9 137 for(Int_t i=0;i<3;i++){
138 fWeighPosSum[i] = 0.;
139 fWeighSigSum[i] = 0.;
140 fAverPosSum[i] = 0.;
141 fWeighPos[i] = 0.;
142 fWeighSig[i] = 0.;
143 fAverPos[i] = 0.;
144 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
145 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
146 fNoEventsContr=0;
147 fTotContributors = 0.;
199775b9 148 }
149}
150//______________________________________________________________________
37add1d1 151Bool_t AliITSMeanVertexer::Init() {
152 // Initialize filters
a3a3a28f 153 // Initialize geometry
37add1d1 154 // Initialize ITS classes
a3a3a28f 155
37add1d1 156 AliGeomManager::LoadGeometry();
157 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f 158
159 AliITSInitGeometry initgeom;
160 AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1 161 if (!geom) return kFALSE;
a3a3a28f 162 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1 163
cf0b1cea 164 fDetTypeRec = new AliITSDetTypeRec();
762133f4 165 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
37add1d1 166 fDetTypeRec->SetITSgeom(geom);
167 fDetTypeRec->SetDefaults();
168 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
169
a3a3a28f 170 // Initialize filter values to their defaults
171 SetFilterOnContributors();
199775b9 172 SetCutOnErrX();
173 SetCutOnR();
174 SetCutOnCls();
a3a3a28f 175
d5eae813 176 // Instatiate vertexer
177 if (!fMode) {
178 fVertexer = new AliITSVertexer3DTapan(1000);
179 }
180 else {
181 fVertexer = new AliITSVertexer3D();
182 fVertexer->SetDetTypeRec(fDetTypeRec);
199775b9 183 AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
4f40f207 184 alias->SetWideFiducialRegion(fZCutDiamond,0.5);
199775b9 185 alias->SetNarrowFiducialRegion(0.5,0.5);
186 alias->SetDeltaPhiCuts(0.5,0.025);
187 alias->SetDCACut(0.1);
188 alias->SetPileupAlgo(3);
189 fVertexer->SetComputeMultiplicity(kFALSE);
d5eae813 190 }
37add1d1 191 return kTRUE;
a3a3a28f 192}
193
194//______________________________________________________________________
195AliITSMeanVertexer::~AliITSMeanVertexer() {
196 // Destructor
37add1d1 197 delete fDetTypeRec;
198 delete fVertexXY;
199 delete fVertexZ;
8e2abaa0 200 delete fMultH;
201 delete fErrXH;
202 delete fMultHa;
203 delete fErrXHa;
204 delete fDistH;
a1b1ad51 205 delete fContrH;
206 delete fContrHa;
d5eae813 207 delete fVertexer;
a3a3a28f 208}
209
210//______________________________________________________________________
d5eae813 211Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
a3a3a28f 212 // Performs SPD local reconstruction
37add1d1 213 // and vertex finding
214 // returns true in case a vertex is found
215
216 // Run SPD cluster finder
8e2abaa0 217 static Int_t evcount = -1;
218 if(evcount <0){
219 evcount++;
220 AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
221 AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
222 }
51a579a9 223// AliCodeTimerAuto("",0);
d5eae813 224 AliITSRecPointContainer::Instance()->PrepareToRead();
37add1d1 225 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
226 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
227
228 Bool_t vtxOK = kFALSE;
199775b9 229 UInt_t nl1=0;
d5eae813 230 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
231 if (!fMode) {
37add1d1 232 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
a3a3a28f 233 }
37add1d1 234 else {
199775b9 235 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
236 nl1=rpcont->GetNClustersInLayerFast(1);
237 if(Filter(vtx,nl1)) vtxOK = kTRUE;
8e2abaa0 238 /* if(vtx){
239 if(vtxOK){
240 printf("The vertex is OK\n");
241 }
242 else {
243 printf("The vertex is NOT OK\n");
244 }
245 vtx->PrintStatus();
246 }
247 else {
248 printf("The vertex was not reconstructed\n");
249 } */
a3a3a28f 250 }
37add1d1 251 delete clustersTree;
199775b9 252 if (vtxOK && fMode){
253 new(fVertArray[fIndex])AliESDVertex(*vtx);
254 fClu0[fIndex]=nl1;
8e2abaa0 255 fAccEvents.SetBitNumber(fIndex);
199775b9 256 fIndex++;
8e2abaa0 257 if(fIndex>=fgkMaxNumOfEvents){
258 if(ComputeMean(kFALSE)){
259 if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
260 }
261 }
199775b9 262 }
37add1d1 263 if (vtx) delete vtx;
a3a3a28f 264
37add1d1 265 return vtxOK;
a3a3a28f 266}
267
268//______________________________________________________________________
37add1d1 269void AliITSMeanVertexer::WriteVertices(const char *filename){
270 // Compute mean vertex and
271 // store it along with the histograms
272 // in a file
273
274 TFile fmv(filename,"update");
199775b9 275 Bool_t itisOK = kFALSE;
276 if(ComputeMean(kFALSE)){
277 if(ComputeMean(kTRUE)){
278 Double_t cov[6];
279 cov[0] = fAverPosSq[0][0]; // variance x
280 cov[1] = fAverPosSq[0][1]; // cov xy
281 cov[2] = fAverPosSq[1][1]; // variance y
282 cov[3] = fAverPosSq[0][2]; // cov xz
283 cov[4] = fAverPosSq[1][2]; // cov yz
284 cov[5] = fAverPosSq[2][2]; // variance z
285 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
286 mv.SetTitle("Mean Vertex");
287 mv.SetName("MeanVertex");
288 // we have to add chi2 here
289 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
37add1d1 290
199775b9 291 mv.Write(mv.GetName(),TObject::kOverwrite);
292 vtx.Write(vtx.GetName(),TObject::kOverwrite);
293 itisOK = kTRUE;
294 }
a3a3a28f 295 }
199775b9 296
297 if(!itisOK){
a3a3a28f 298 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
299 }
37add1d1 300
301 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
302 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
8e2abaa0 303 fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
304 fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
305 fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
306 fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
307 fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
a1b1ad51 308 fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
309 fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
37add1d1 310 fmv.Close();
a3a3a28f 311}
312
313//______________________________________________________________________
199775b9 314Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
a3a3a28f 315 // Apply selection criteria to events
316 Bool_t status = kFALSE;
199775b9 317 if(!vert)return status;
8d7f6205 318 // Remove vertices reconstructed with vertexerZ
319 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
a3a3a28f 320 Int_t ncontr = vert->GetNContributors();
a3a3a28f 321 AliDebug(1,Form("Number of contributors = %d",ncontr));
199775b9 322 Double_t ex = vert->GetXRes();
8e2abaa0 323 fMultH->Fill(mult);
324 fErrXH->Fill(ex*1000.);
a1b1ad51 325 fContrH->Fill((Float_t)ncontr);
199775b9 326 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
327 ((ex*1000.)<fErrXCut)) {
328 status = kTRUE;
8e2abaa0 329 fMultHa->Fill(mult);
330 fErrXHa->Fill(ex*1000.);
a1b1ad51 331 fContrHa->Fill((Float_t)ncontr);
199775b9 332 }
a3a3a28f 333 return status;
334}
335
336//______________________________________________________________________
337void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
338 // update mean vertex
339 Double_t currentPos[3],currentSigma[3];
340 vert->GetXYZ(currentPos);
341 vert->GetSigmaXYZ(currentSigma);
342 Bool_t goon = kTRUE;
343 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
344 if(!goon)return;
345 for(Int_t i=0;i<3;i++){
762133f4 346 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
347 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
348 fAverPosSum[i]+=currentPos[i];
a3a3a28f 349 }
350 for(Int_t i=0;i<3;i++){
351 for(Int_t j=i;j<3;j++){
762133f4 352 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f 353 }
354 }
37add1d1 355
199775b9 356 fTotContributors+=vert->GetNContributors();
37add1d1 357 fVertexXY->Fill(currentPos[0],currentPos[1]);
358 fVertexZ->Fill(currentPos[2]);
359
a3a3a28f 360 fNoEventsContr++;
361}
362
363//______________________________________________________________________
199775b9 364 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
365 // compute mean vertex
366 // called once with killOutliers = kFALSE and then again with
367 // killOutliers = kTRUE
8e2abaa0 368
a1b1ad51 369 static Bool_t preliminary = kTRUE;
370 static Int_t oldnumbevt = 0;
371 if(!(preliminary || killOutliers))return kTRUE; //ComputeMean is
8e2abaa0 372 // executed only once with argument kFALSE
199775b9 373 Double_t wpos[3];
374 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
8e2abaa0 375
a1b1ad51 376 Int_t istart = oldnumbevt;
377 if(killOutliers)istart = 0;
378 if(killOutliers && preliminary){
379 preliminary = kFALSE;
380 Reset(kTRUE,kFALSE);
8e2abaa0 381 }
a1b1ad51 382 oldnumbevt = fVertArray.GetEntries();
383 for(Int_t i=istart;i<fVertArray.GetEntries();i++){
199775b9 384 AliESDVertex* vert = NULL;
762133f4 385
199775b9 386 // second pass
387 if(killOutliers && fAccEvents.TestBitNumber(i)){
388 vert=(AliESDVertex*)fVertArray[i];
389 Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
390 dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
391 dist=sqrt(dist)*10.; // distance in mm
8e2abaa0 392 fDistH->Fill(dist);
199775b9 393 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
394 }
395
396 if(!fAccEvents.TestBitNumber(i))continue;
397 vert=(AliESDVertex*)fVertArray[i];
398 AddToMean(vert);
399 }
a1b1ad51 400 Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
401 if(bad) {
402 if(killOutliers){
403// when the control reachs this point, the preliminary evaluation of the
404// diamond region has to be redone. So a general reset must be issued
405// if this occurs, it is likely that the cuts are badly defined
406 ResetArray();
407 Reset(kFALSE,kTRUE);
408 preliminary = kTRUE;
409 oldnumbevt = 0;
410 }
411 return kFALSE;
412 }
199775b9 413 Double_t nevents = fNoEventsContr;
414 for(Int_t i=0;i<3;i++){
415 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
416 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
417 fAverPos[i] = fAverPosSum[i]/nevents;
418 }
419 for(Int_t i=0;i<3;i++){
420 for(Int_t j=i;j<3;j++){
421 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
422 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
423 }
424 }
8e2abaa0 425 if(killOutliers)ResetArray();
199775b9 426 fAverContributors = fTotContributors/nevents;
427 return kTRUE;
428 }