5 #include "AliGeomManager.h"
6 #include "AliITSDetTypeRec.h"
7 #include "AliITSInitGeometry.h"
8 #include "AliITSMeanVertexer.h"
9 #include "AliITSRecPointContainer.h"
10 #include "AliITSLoader.h"
12 #include "AliRawReader.h"
13 #include "AliRawReaderDate.h"
14 #include "AliRawReaderRoot.h"
15 #include "AliRunLoader.h"
16 #include "AliITSVertexer3D.h"
17 #include "AliITSVertexer3DTapan.h"
18 #include "AliESDVertex.h"
19 #include "AliMeanVertex.h"
20 //#include "AliCodeTimer.h"
22 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
23 ClassImp(AliITSMeanVertexer)
24 ///////////////////////////////////////////////////////////////////////
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) //
32 // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
33 ///////////////////////////////////////////////////////////////////////
35 //______________________________________________________________________
36 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
42 fAverContributors(0.),
43 fFilterOnContributors(0),
46 fAccEvents(fgkMaxNumOfEvents),
47 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
60 // Default Constructor
63 // Histograms initialization
64 const Float_t xLimit = 0.6, yLimit = 0.6, zLimit = 50.0;
65 const Float_t xDelta = 0.003, yDelta = 0.003, zDelta = 0.2;
66 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
67 2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
68 2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
69 fVertexXY->SetXTitle("X , cm");
70 fVertexXY->SetYTitle("Y , cm");
71 fVertexXY->SetOption("colz");
72 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
73 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
74 fVertexZ->SetXTitle("Z , cm");
75 fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
76 fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
77 fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
78 fErrXHa->SetLineColor(kRed);
79 fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
80 fMultHa->SetLineColor(kRed);
81 fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
82 fClu0 = new UInt_t [fgkMaxNumOfEvents];
83 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
84 fAccEvents.ResetAllBits(kTRUE);
87 //______________________________________________________________________
88 void AliITSMeanVertexer::Reset(){
90 static Int_t counter = 0;
93 Double_t rangeX=fVertexXY->GetXaxis()->GetXmax()-fVertexXY->GetXaxis()->GetXmin();
94 Double_t rangeY=fVertexXY->GetYaxis()->GetXmax()-fVertexXY->GetYaxis()->GetXmin();
95 Int_t nx=fVertexXY->GetNbinsX();
96 Int_t ny=fVertexXY->GetNbinsY();
98 Double_t xmi=fWeighPos[0]-rangeX/2.;
99 Double_t xma=fWeighPos[0]+rangeX/2.;
100 Double_t ymi=fWeighPos[1]-rangeY/2.;
101 Double_t yma=fWeighPos[1]+rangeY/2.;
102 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",nx,xmi,xma,ny,ymi,yma);
103 fVertexXY->SetXTitle("X , cm");
104 fVertexXY->SetYTitle("Y , cm");
105 fVertexXY->SetOption("colz");
106 fVertexXY->SetDirectory(0);
108 for(Int_t i=0;i<3;i++){
109 fWeighPosSum[i] = 0.;
110 fWeighSigSum[i] = 0.;
115 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
116 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
118 fTotContributors = 0.;
119 if(fVertexZ)fVertexZ->Reset();
122 //______________________________________________________________________
123 Bool_t AliITSMeanVertexer::Init() {
124 // Initialize filters
125 // Initialize geometry
126 // Initialize ITS classes
128 AliGeomManager::LoadGeometry();
129 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
131 AliITSInitGeometry initgeom;
132 AliITSgeom *geom = initgeom.CreateAliITSgeom();
133 if (!geom) return kFALSE;
134 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
136 fDetTypeRec = new AliITSDetTypeRec();
137 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
138 fDetTypeRec->SetITSgeom(geom);
139 fDetTypeRec->SetDefaults();
140 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
142 // Initialize filter values to their defaults
143 SetFilterOnContributors();
148 // Instatiate vertexer
150 fVertexer = new AliITSVertexer3DTapan(1000);
153 fVertexer = new AliITSVertexer3D();
154 fVertexer->SetDetTypeRec(fDetTypeRec);
155 AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
156 alias->SetWideFiducialRegion(40.,0.5);
157 alias->SetNarrowFiducialRegion(0.5,0.5);
158 alias->SetDeltaPhiCuts(0.5,0.025);
159 alias->SetDCACut(0.1);
160 alias->SetPileupAlgo(3);
161 fVertexer->SetComputeMultiplicity(kFALSE);
166 //______________________________________________________________________
167 AliITSMeanVertexer::~AliITSMeanVertexer() {
180 //______________________________________________________________________
181 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
182 // Performs SPD local reconstruction
183 // and vertex finding
184 // returns true in case a vertex is found
186 // Run SPD cluster finder
187 static Int_t evcount = -1;
190 AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
191 AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
193 // AliCodeTimerAuto("",0);
194 AliITSRecPointContainer::Instance()->PrepareToRead();
195 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
196 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
198 Bool_t vtxOK = kFALSE;
200 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
202 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
205 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
206 nl1=rpcont->GetNClustersInLayerFast(1);
207 if(Filter(vtx,nl1)) vtxOK = kTRUE;
210 printf("The vertex is OK\n");
213 printf("The vertex is NOT OK\n");
218 printf("The vertex was not reconstructed\n");
223 new(fVertArray[fIndex])AliESDVertex(*vtx);
225 fAccEvents.SetBitNumber(fIndex);
227 if(fIndex>=fgkMaxNumOfEvents){
228 if(ComputeMean(kFALSE)){
229 if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
238 //______________________________________________________________________
239 void AliITSMeanVertexer::WriteVertices(const char *filename){
240 // Compute mean vertex and
241 // store it along with the histograms
244 TFile fmv(filename,"update");
245 Bool_t itisOK = kFALSE;
246 if(ComputeMean(kFALSE)){
247 if(ComputeMean(kTRUE)){
249 cov[0] = fAverPosSq[0][0]; // variance x
250 cov[1] = fAverPosSq[0][1]; // cov xy
251 cov[2] = fAverPosSq[1][1]; // variance y
252 cov[3] = fAverPosSq[0][2]; // cov xz
253 cov[4] = fAverPosSq[1][2]; // cov yz
254 cov[5] = fAverPosSq[2][2]; // variance z
255 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
256 mv.SetTitle("Mean Vertex");
257 mv.SetName("MeanVertex");
258 // we have to add chi2 here
259 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
261 mv.Write(mv.GetName(),TObject::kOverwrite);
262 vtx.Write(vtx.GetName(),TObject::kOverwrite);
268 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
271 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
272 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
273 fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
274 fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
275 fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
276 fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
277 fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
281 //______________________________________________________________________
282 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
283 // Apply selection criteria to events
284 Bool_t status = kFALSE;
285 if(!vert)return status;
286 // Remove vertices reconstructed with vertexerZ
287 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
288 Int_t ncontr = vert->GetNContributors();
289 AliDebug(1,Form("Number of contributors = %d",ncontr));
290 Double_t ex = vert->GetXRes();
292 fErrXH->Fill(ex*1000.);
293 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
294 ((ex*1000.)<fErrXCut)) {
297 fErrXHa->Fill(ex*1000.);
302 //______________________________________________________________________
303 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
304 // update mean vertex
305 Double_t currentPos[3],currentSigma[3];
306 vert->GetXYZ(currentPos);
307 vert->GetSigmaXYZ(currentSigma);
309 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
311 for(Int_t i=0;i<3;i++){
312 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
313 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
314 fAverPosSum[i]+=currentPos[i];
316 for(Int_t i=0;i<3;i++){
317 for(Int_t j=i;j<3;j++){
318 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
322 fTotContributors+=vert->GetNContributors();
323 fVertexXY->Fill(currentPos[0],currentPos[1]);
324 fVertexZ->Fill(currentPos[2]);
329 //______________________________________________________________________
330 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
331 // compute mean vertex
332 // called once with killOutliers = kFALSE and then again with
333 // killOutliers = kTRUE
335 static Bool_t firstime = kTRUE;
336 if(fNoEventsContr>2 && !killOutliers)return kTRUE; //ComputeMean is
337 // executed only once with argument kFALSE
339 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
341 if(killOutliers && firstime){
346 for(Int_t i=0;i<fVertArray.GetEntries();i++){
347 AliESDVertex* vert = NULL;
350 if(killOutliers && fAccEvents.TestBitNumber(i)){
351 vert=(AliESDVertex*)fVertArray[i];
352 Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
353 dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
354 dist=sqrt(dist)*10.; // distance in mm
356 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
359 if(!fAccEvents.TestBitNumber(i))continue;
360 vert=(AliESDVertex*)fVertArray[i];
363 if(fNoEventsContr < 2) return kFALSE;
364 Double_t nevents = fNoEventsContr;
365 for(Int_t i=0;i<3;i++){
366 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
367 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
368 fAverPos[i] = fAverPosSum[i]/nevents;
370 for(Int_t i=0;i<3;i++){
371 for(Int_t j=i;j<3;j++){
372 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
373 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
376 if(killOutliers)ResetArray();
377 fAverContributors = fTotContributors/nevents;