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)
25 ///////////////////////////////////////////////////////////////////////
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) //
33 // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
34 ///////////////////////////////////////////////////////////////////////
38 //______________________________________________________________________
39 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
45 fAverContributors(0.),
46 fFilterOnContributors(0),
49 fAccEvents(fgkMaxNumOfEvents),
50 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
63 // Default Constructor
66 // Histograms initialization
67 const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
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);
72 fVertexXY->SetXTitle("X , cm");
73 fVertexXY->SetYTitle("Y , cm");
74 fVertexXY->SetOption("colz");
75 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
76 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
77 fVertexZ->SetXTitle("Z , cm");
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.);
81 fErrXHa->SetLineColor(kRed);
82 fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
83 fMultHa->SetLineColor(kRed);
84 fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
85 fClu0 = new UInt_t [fgkMaxNumOfEvents];
86 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
87 fAccEvents.ResetAllBits(kTRUE);
90 //______________________________________________________________________
91 void AliITSMeanVertexer::Reset(){
93 for(Int_t i=0;i<3;i++){
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.;
103 fTotContributors = 0.;
104 if(fVertexXY)fVertexXY->Reset();
105 if(fVertexZ)fVertexZ->Reset();
108 //______________________________________________________________________
109 Bool_t AliITSMeanVertexer::Init() {
110 // Initialize filters
111 // Initialize geometry
112 // Initialize ITS classes
114 AliGeomManager::LoadGeometry();
115 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
117 AliITSInitGeometry initgeom;
118 AliITSgeom *geom = initgeom.CreateAliITSgeom();
119 if (!geom) return kFALSE;
120 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
122 fDetTypeRec = new AliITSDetTypeRec();
123 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
124 fDetTypeRec->SetITSgeom(geom);
125 fDetTypeRec->SetDefaults();
126 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
128 // Initialize filter values to their defaults
129 SetFilterOnContributors();
134 // Instatiate vertexer
136 fVertexer = new AliITSVertexer3DTapan(1000);
139 fVertexer = new AliITSVertexer3D();
140 fVertexer->SetDetTypeRec(fDetTypeRec);
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);
152 //______________________________________________________________________
153 AliITSMeanVertexer::~AliITSMeanVertexer() {
166 //______________________________________________________________________
167 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
168 // Performs SPD local reconstruction
169 // and vertex finding
170 // returns true in case a vertex is found
172 // Run SPD cluster finder
173 static Int_t evcount = -1;
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));
179 // AliCodeTimerAuto("",0);
180 AliITSRecPointContainer::Instance()->PrepareToRead();
181 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
182 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
184 Bool_t vtxOK = kFALSE;
186 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
188 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
191 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
192 nl1=rpcont->GetNClustersInLayerFast(1);
193 if(Filter(vtx,nl1)) vtxOK = kTRUE;
196 printf("The vertex is OK\n");
199 printf("The vertex is NOT OK\n");
204 printf("The vertex was not reconstructed\n");
209 new(fVertArray[fIndex])AliESDVertex(*vtx);
211 fAccEvents.SetBitNumber(fIndex);
213 if(fIndex>=fgkMaxNumOfEvents){
214 if(ComputeMean(kFALSE)){
215 if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
224 //______________________________________________________________________
225 void AliITSMeanVertexer::WriteVertices(const char *filename){
226 // Compute mean vertex and
227 // store it along with the histograms
230 TFile fmv(filename,"update");
231 Bool_t itisOK = kFALSE;
232 if(ComputeMean(kFALSE)){
233 if(ComputeMean(kTRUE)){
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");
247 mv.Write(mv.GetName(),TObject::kOverwrite);
248 vtx.Write(vtx.GetName(),TObject::kOverwrite);
254 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
257 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
258 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
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);
267 //______________________________________________________________________
268 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
269 // Apply selection criteria to events
270 Bool_t status = kFALSE;
271 if(!vert)return status;
272 // Remove vertices reconstructed with vertexerZ
273 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
274 Int_t ncontr = vert->GetNContributors();
275 AliDebug(1,Form("Number of contributors = %d",ncontr));
276 Double_t ex = vert->GetXRes();
278 fErrXH->Fill(ex*1000.);
279 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
280 ((ex*1000.)<fErrXCut)) {
283 fErrXHa->Fill(ex*1000.);
288 //______________________________________________________________________
289 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
290 // update mean vertex
291 Double_t currentPos[3],currentSigma[3];
292 vert->GetXYZ(currentPos);
293 vert->GetSigmaXYZ(currentSigma);
295 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
297 for(Int_t i=0;i<3;i++){
298 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
299 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
300 fAverPosSum[i]+=currentPos[i];
302 for(Int_t i=0;i<3;i++){
303 for(Int_t j=i;j<3;j++){
304 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
308 fTotContributors+=vert->GetNContributors();
309 fVertexXY->Fill(currentPos[0],currentPos[1]);
310 fVertexZ->Fill(currentPos[2]);
315 //______________________________________________________________________
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
321 static Bool_t firstime = kTRUE;
322 if(fNoEventsContr>2 && !killOutliers)return kTRUE; //ComputeMean is
323 // executed only once with argument kFALSE
325 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
327 if(killOutliers && firstime){
332 for(Int_t i=0;i<fVertArray.GetEntries();i++){
333 AliESDVertex* vert = NULL;
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
342 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
345 if(!fAccEvents.TestBitNumber(i))continue;
346 vert=(AliESDVertex*)fVertArray[i];
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;
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];
362 if(killOutliers)ResetArray();
363 fAverContributors = fTotContributors/nevents;