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