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