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),
63 // Default Constructor
64 SetZFiducialRegion(); // sets fZCutDiamond to its default value
67 // Histograms initialization
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;
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);
73 fVertexXY->SetXTitle("X , cm");
74 fVertexXY->SetYTitle("Y , cm");
75 fVertexXY->SetOption("colz");
76 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
77 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
78 fVertexZ->SetXTitle("Z , cm");
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.);
82 fErrXHa->SetLineColor(kRed);
83 fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
84 fMultHa->SetLineColor(kRed);
85 fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
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);
90 fClu0 = new UInt_t [fgkMaxNumOfEvents];
91 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
92 fAccEvents.ResetAllBits(kTRUE);
95 //______________________________________________________________________
96 void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
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
101 // redefine2D is FALSE and complete is TRUE when the method is called
102 // after a failure cof ComputeMean(kTRUE)
104 static Int_t counter = 0;
105 if(fVertexXY && redefine2D){
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();
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);
122 else if(fVertexXY && !redefine2D){
137 for(Int_t i=0;i<3;i++){
138 fWeighPosSum[i] = 0.;
139 fWeighSigSum[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.;
147 fTotContributors = 0.;
150 //______________________________________________________________________
151 Bool_t AliITSMeanVertexer::Init() {
152 // Initialize filters
153 // Initialize geometry
154 // Initialize ITS classes
156 AliGeomManager::LoadGeometry();
157 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
159 AliITSInitGeometry initgeom;
160 AliITSgeom *geom = initgeom.CreateAliITSgeom();
161 if (!geom) return kFALSE;
162 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
164 fDetTypeRec = new AliITSDetTypeRec();
165 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
166 fDetTypeRec->SetITSgeom(geom);
167 fDetTypeRec->SetDefaults();
168 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
170 // Initialize filter values to their defaults
171 SetFilterOnContributors();
176 // Instatiate vertexer
178 fVertexer = new AliITSVertexer3DTapan(1000);
181 fVertexer = new AliITSVertexer3D();
182 fVertexer->SetDetTypeRec(fDetTypeRec);
183 AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
184 alias->SetWideFiducialRegion(fZCutDiamond,0.5);
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);
194 //______________________________________________________________________
195 AliITSMeanVertexer::~AliITSMeanVertexer() {
210 //______________________________________________________________________
211 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
212 // Performs SPD local reconstruction
213 // and vertex finding
214 // returns true in case a vertex is found
216 // Run SPD cluster finder
217 static Int_t evcount = -1;
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));
223 // AliCodeTimerAuto("",0);
224 AliITSRecPointContainer::Instance()->PrepareToRead();
225 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
226 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
228 Bool_t vtxOK = kFALSE;
230 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
232 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
235 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
236 nl1=rpcont->GetNClustersInLayerFast(1);
237 if(Filter(vtx,nl1)) vtxOK = kTRUE;
240 printf("The vertex is OK\n");
243 printf("The vertex is NOT OK\n");
248 printf("The vertex was not reconstructed\n");
253 new(fVertArray[fIndex])AliESDVertex(*vtx);
255 fAccEvents.SetBitNumber(fIndex);
257 if(fIndex>=fgkMaxNumOfEvents){
258 if(ComputeMean(kFALSE)){
259 if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
268 //______________________________________________________________________
269 void AliITSMeanVertexer::WriteVertices(const char *filename){
270 // Compute mean vertex and
271 // store it along with the histograms
274 TFile fmv(filename,"update");
275 Bool_t itisOK = kFALSE;
276 if(ComputeMean(kFALSE)){
277 if(ComputeMean(kTRUE)){
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 // We use standard average and not weighed averag now
286 // AliMeanVertex is apparently not taken by the preprocessor; only
287 // the AliESDVertex object is retrieved
288 // AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
289 AliMeanVertex mv(fAverPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
290 mv.SetTitle("Mean Vertex");
291 mv.SetName("MeanVertex");
292 // we have to add chi2 here
293 // AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
294 AliESDVertex vtx(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
296 mv.Write(mv.GetName(),TObject::kOverwrite);
297 vtx.Write(vtx.GetName(),TObject::kOverwrite);
303 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
306 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
307 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
308 fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
309 fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
310 fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
311 fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
312 fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
313 fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
314 fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
318 //______________________________________________________________________
319 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
320 // Apply selection criteria to events
321 Bool_t status = kFALSE;
322 if(!vert)return status;
323 // Remove vertices reconstructed with vertexerZ
324 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
325 Int_t ncontr = vert->GetNContributors();
326 AliDebug(1,Form("Number of contributors = %d",ncontr));
327 Double_t ex = vert->GetXRes();
329 fErrXH->Fill(ex*1000.);
330 fContrH->Fill((Float_t)ncontr);
331 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
332 ((ex*1000.)<fErrXCut)) {
335 fErrXHa->Fill(ex*1000.);
336 fContrHa->Fill((Float_t)ncontr);
341 //______________________________________________________________________
342 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
343 // update mean vertex
344 Double_t currentPos[3],currentSigma[3];
345 vert->GetXYZ(currentPos);
346 vert->GetSigmaXYZ(currentSigma);
348 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
350 for(Int_t i=0;i<3;i++){
351 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
352 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
353 fAverPosSum[i]+=currentPos[i];
355 for(Int_t i=0;i<3;i++){
356 for(Int_t j=i;j<3;j++){
357 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
361 fTotContributors+=vert->GetNContributors();
362 fVertexXY->Fill(currentPos[0],currentPos[1]);
363 fVertexZ->Fill(currentPos[2]);
368 //______________________________________________________________________
369 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
370 // compute mean vertex
371 // called once with killOutliers = kFALSE and then again with
372 // killOutliers = kTRUE
374 static Bool_t preliminary = kTRUE;
375 static Int_t oldnumbevt = 0;
376 if(!(preliminary || killOutliers))return kTRUE; //ComputeMean is
377 // executed only once with argument kFALSE
379 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
381 Int_t istart = oldnumbevt;
382 if(killOutliers)istart = 0;
383 if(killOutliers && preliminary){
384 preliminary = kFALSE;
387 oldnumbevt = fVertArray.GetEntries();
388 for(Int_t i=istart;i<fVertArray.GetEntries();i++){
389 AliESDVertex* vert = NULL;
392 if(killOutliers && fAccEvents.TestBitNumber(i)){
393 vert=(AliESDVertex*)fVertArray[i];
394 Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
395 dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
396 dist=sqrt(dist)*10.; // distance in mm
398 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
401 if(!fAccEvents.TestBitNumber(i))continue;
402 vert=(AliESDVertex*)fVertArray[i];
405 Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
408 // when the control reachs this point, the preliminary evaluation of the
409 // diamond region has to be redone. So a general reset must be issued
410 // if this occurs, it is likely that the cuts are badly defined
418 Double_t nevents = fNoEventsContr;
419 for(Int_t i=0;i<3;i++){
420 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
421 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
422 fAverPos[i] = fAverPosSum[i]/nevents;
424 for(Int_t i=0;i<3;i++){
425 for(Int_t j=i;j<3;j++){
426 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
427 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
430 if(killOutliers)ResetArray();
431 fAverContributors = fTotContributors/nevents;