Line data Source code
1 : #include <TFile.h>
2 : #include <TH1.h>
3 : #include <TH2.h>
4 : #include <TTree.h>
5 : #include "AliGeomManager.h"
6 : #include "AliITSDetTypeRec.h"
7 : #include "AliITSInitGeometry.h"
8 : #include "AliITSMeanVertexer.h"
9 : #include "AliITSRecPointContainer.h"
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"
17 : #include "AliITSVertexer3DTapan.h"
18 : #include "AliESDVertex.h"
19 : #include "AliMeanVertex.h"
20 : //#include "AliCodeTimer.h"
21 :
22 : const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
23 116 : ClassImp(AliITSMeanVertexer)
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 : // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
33 : ///////////////////////////////////////////////////////////////////////
34 : /* $Id$ */
35 : //______________________________________________________________________
36 0 : AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
37 0 : fDetTypeRec(NULL),
38 0 : fVertexXY(NULL),
39 0 : fVertexZ(NULL),
40 0 : fNoEventsContr(0),
41 0 : fTotContributors(0.),
42 0 : fAverContributors(0.),
43 0 : fFilterOnContributors(0),
44 0 : fMode(mode),
45 0 : fVertexer(NULL),
46 0 : fAccEvents(fgkMaxNumOfEvents),
47 0 : fVertArray("AliESDVertex",fgkMaxNumOfEvents),
48 0 : fClu0(NULL),
49 0 : fIndex(0),
50 0 : fErrXCut(0.),
51 0 : fRCut(0.),
52 0 : fZCutDiamond(0.),
53 0 : fLowSPD0(0),
54 0 : fHighSPD0(0),
55 0 : fMultH(NULL),
56 0 : fErrXH(NULL),
57 0 : fMultHa(NULL),
58 0 : fErrXHa(NULL),
59 0 : fDistH(NULL),
60 0 : fContrH(NULL),
61 0 : fContrHa(NULL)
62 0 : {
63 : // Default Constructor
64 0 : SetZFiducialRegion(); // sets fZCutDiamond to its default value
65 0 : Reset(kFALSE,kFALSE);
66 :
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 0 : 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 0 : fVertexXY->SetXTitle("X , cm");
74 0 : fVertexXY->SetYTitle("Y , cm");
75 0 : fVertexXY->SetOption("colz");
76 0 : fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
77 : 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
78 0 : fVertexZ->SetXTitle("Z , cm");
79 0 : fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
80 0 : fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
81 0 : fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
82 0 : fErrXHa->SetLineColor(kRed);
83 0 : fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
84 0 : fMultHa->SetLineColor(kRed);
85 0 : fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
86 0 : fContrH = new TH1F("contrib","Number of contributors - before cuts",300,0.5,300.5);
87 0 : fContrHa = new TH1F("contriba","Number of contributors - after cuts",300,0.5,300.5);
88 0 : fContrHa->SetLineColor(kRed);
89 :
90 0 : fClu0 = new UInt_t [fgkMaxNumOfEvents];
91 0 : for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
92 0 : fAccEvents.ResetAllBits(kTRUE);
93 0 : }
94 :
95 : //______________________________________________________________________
96 : void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
97 : // Reset averages
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
100 : // pass.
101 : // redefine2D is FALSE and complete is TRUE when the method is called
102 : // after a failure cof ComputeMean(kTRUE)
103 :
104 : static Int_t counter = 0;
105 0 : if(fVertexXY && redefine2D){
106 0 : counter++;
107 0 : Double_t rangeX=fVertexXY->GetXaxis()->GetXmax()-fVertexXY->GetXaxis()->GetXmin();
108 0 : Double_t rangeY=fVertexXY->GetYaxis()->GetXmax()-fVertexXY->GetYaxis()->GetXmin();
109 0 : Int_t nx=fVertexXY->GetNbinsX();
110 0 : Int_t ny=fVertexXY->GetNbinsY();
111 0 : delete fVertexXY;
112 0 : Double_t xmi=fWeighPos[0]-rangeX/2.;
113 0 : Double_t xma=fWeighPos[0]+rangeX/2.;
114 0 : Double_t ymi=fWeighPos[1]-rangeY/2.;
115 0 : Double_t yma=fWeighPos[1]+rangeY/2.;
116 0 : fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",nx,xmi,xma,ny,ymi,yma);
117 0 : fVertexXY->SetXTitle("X , cm");
118 0 : fVertexXY->SetYTitle("Y , cm");
119 0 : fVertexXY->SetOption("colz");
120 0 : fVertexXY->SetDirectory(0);
121 0 : }
122 0 : else if(fVertexXY && !redefine2D){
123 0 : fVertexXY->Reset();
124 0 : }
125 0 : if(fVertexZ){
126 0 : fVertexZ->Reset();
127 0 : fDistH->Reset();
128 0 : if(complete){
129 0 : fErrXH->Reset();
130 0 : fMultH->Reset();
131 0 : fErrXHa->Reset();
132 0 : fMultHa->Reset();
133 0 : fContrH->Reset();
134 0 : fContrHa->Reset();
135 0 : }
136 : }
137 0 : for(Int_t i=0;i<3;i++){
138 0 : fWeighPosSum[i] = 0.;
139 0 : fWeighSigSum[i] = 0.;
140 0 : fAverPosSum[i] = 0.;
141 0 : fWeighPos[i] = 0.;
142 0 : fWeighSig[i] = 0.;
143 0 : fAverPos[i] = 0.;
144 0 : for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
145 0 : for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
146 0 : fNoEventsContr=0;
147 0 : fTotContributors = 0.;
148 : }
149 0 : }
150 : //______________________________________________________________________
151 : Bool_t AliITSMeanVertexer::Init() {
152 : // Initialize filters
153 : // Initialize geometry
154 : // Initialize ITS classes
155 :
156 0 : AliGeomManager::LoadGeometry();
157 0 : if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
158 :
159 0 : AliITSInitGeometry initgeom;
160 0 : AliITSgeom *geom = initgeom.CreateAliITSgeom();
161 0 : if (!geom) return kFALSE;
162 0 : printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
163 :
164 0 : fDetTypeRec = new AliITSDetTypeRec();
165 0 : fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
166 0 : fDetTypeRec->SetITSgeom(geom);
167 0 : fDetTypeRec->SetDefaults();
168 0 : fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
169 :
170 : // Initialize filter values to their defaults
171 0 : SetFilterOnContributors();
172 0 : SetCutOnErrX();
173 0 : SetCutOnR();
174 0 : SetCutOnCls();
175 :
176 : // Instatiate vertexer
177 0 : if (!fMode) {
178 0 : fVertexer = new AliITSVertexer3DTapan(1000);
179 0 : }
180 : else {
181 0 : fVertexer = new AliITSVertexer3D();
182 0 : fVertexer->SetDetTypeRec(fDetTypeRec);
183 0 : AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
184 : // alias->SetWideFiducialRegion(fZCutDiamond,0.5);
185 0 : alias->SetWideFiducialRegion(fZCutDiamond,2.5);
186 0 : alias->SetNarrowFiducialRegion(0.5,0.5);
187 0 : alias->SetDeltaPhiCuts(0.5,0.025);
188 0 : alias->SetDCACut(0.1);
189 0 : alias->SetPileupAlgo(3);
190 0 : alias->SetFallBack(6000);
191 0 : fVertexer->SetComputeMultiplicity(kFALSE);
192 : }
193 0 : return kTRUE;
194 0 : }
195 :
196 : //______________________________________________________________________
197 0 : AliITSMeanVertexer::~AliITSMeanVertexer() {
198 : // Destructor
199 0 : delete fDetTypeRec;
200 0 : delete fVertexXY;
201 0 : delete fVertexZ;
202 0 : delete fMultH;
203 0 : delete fErrXH;
204 0 : delete fMultHa;
205 0 : delete fErrXHa;
206 0 : delete fDistH;
207 0 : delete fContrH;
208 0 : delete fContrHa;
209 0 : delete fVertexer;
210 0 : }
211 :
212 : //______________________________________________________________________
213 : Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
214 : // Performs SPD local reconstruction
215 : // and vertex finding
216 : // returns true in case a vertex is found
217 :
218 : // Run SPD cluster finder
219 : static Int_t evcount = -1;
220 0 : if(evcount <0){
221 0 : evcount++;
222 0 : AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
223 0 : AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
224 0 : }
225 : // AliCodeTimerAuto("",0);
226 0 : AliITSRecPointContainer::Instance()->PrepareToRead();
227 0 : TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
228 0 : fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
229 :
230 : Bool_t vtxOK = kFALSE;
231 : UInt_t nl1=0;
232 0 : nl1=AliITSRecPointContainer::Instance()->GetNClustersInLayerFast(1);
233 0 : fMultH->Fill(nl1);
234 0 : if(!(nl1>fLowSPD0 && nl1<fHighSPD0)){
235 0 : delete clustersTree;
236 0 : return kFALSE;
237 : }
238 0 : AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
239 0 : if (!fMode) {
240 0 : if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
241 : }
242 : else {
243 0 : AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
244 0 : nl1=rpcont->GetNClustersInLayerFast(1);
245 0 : if(Filter(vtx,nl1)) vtxOK = kTRUE;
246 : /* if(vtx){
247 : if(vtxOK){
248 : printf("The vertex is OK\n");
249 : }
250 : else {
251 : printf("The vertex is NOT OK\n");
252 : }
253 : vtx->PrintStatus();
254 : }
255 : else {
256 : printf("The vertex was not reconstructed\n");
257 : } */
258 : }
259 0 : delete clustersTree;
260 0 : if (vtxOK && fMode){
261 0 : new(fVertArray[fIndex])AliESDVertex(*vtx);
262 0 : fClu0[fIndex]=nl1;
263 0 : fAccEvents.SetBitNumber(fIndex);
264 0 : fIndex++;
265 0 : if(fIndex>=fgkMaxNumOfEvents){
266 0 : if(ComputeMean(kFALSE)){
267 0 : if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
268 : }
269 : }
270 : }
271 0 : if (vtx) delete vtx;
272 :
273 : return vtxOK;
274 0 : }
275 :
276 : //______________________________________________________________________
277 : void AliITSMeanVertexer::WriteVertices(const char *filename){
278 : // Compute mean vertex and
279 : // store it along with the histograms
280 : // in a file
281 :
282 0 : TFile fmv(filename,"update");
283 : Bool_t itisOK = kFALSE;
284 0 : if(ComputeMean(kFALSE)){
285 0 : if(ComputeMean(kTRUE)){
286 0 : Double_t cov[6];
287 0 : cov[0] = fAverPosSq[0][0]; // variance x
288 0 : cov[1] = fAverPosSq[0][1]; // cov xy
289 0 : cov[2] = fAverPosSq[1][1]; // variance y
290 0 : cov[3] = fAverPosSq[0][2]; // cov xz
291 0 : cov[4] = fAverPosSq[1][2]; // cov yz
292 0 : cov[5] = fAverPosSq[2][2]; // variance z
293 : // We use standard average and not weighed averag now
294 : // AliMeanVertex is apparently not taken by the preprocessor; only
295 : // the AliESDVertex object is retrieved
296 : // AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
297 0 : AliMeanVertex mv(fAverPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
298 0 : mv.SetTitle("Mean Vertex");
299 0 : mv.SetName("MeanVertex");
300 : // we have to add chi2 here
301 : // AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
302 0 : AliESDVertex vtx(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
303 :
304 0 : mv.Write(mv.GetName(),TObject::kOverwrite);
305 0 : vtx.Write(vtx.GetName(),TObject::kOverwrite);
306 : itisOK = kTRUE;
307 0 : }
308 : }
309 :
310 0 : if(!itisOK){
311 0 : AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
312 : }
313 :
314 0 : fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
315 0 : fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
316 0 : fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
317 0 : fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
318 0 : fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
319 0 : fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
320 0 : fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
321 0 : fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
322 0 : fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
323 0 : fmv.Close();
324 0 : }
325 :
326 : //______________________________________________________________________
327 : Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
328 : // Apply selection criteria to events
329 : Bool_t status = kFALSE;
330 0 : if(!vert)return status;
331 : // Remove vertices reconstructed with vertexerZ
332 0 : if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
333 0 : Int_t ncontr = vert->GetNContributors();
334 0 : AliDebug(1,Form("Number of contributors = %d",ncontr));
335 0 : Double_t ex = vert->GetXRes();
336 : // fMultH->Fill(mult);
337 0 : fErrXH->Fill(ex*1000.);
338 0 : fContrH->Fill((Float_t)ncontr);
339 0 : if(ncontr>fFilterOnContributors && ((ex*1000.)<fErrXCut)) {
340 : status = kTRUE;
341 0 : fMultHa->Fill(mult);
342 0 : fErrXHa->Fill(ex*1000.);
343 0 : fContrHa->Fill((Float_t)ncontr);
344 0 : }
345 0 : return status;
346 0 : }
347 :
348 : //______________________________________________________________________
349 : void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
350 : // update mean vertex
351 0 : Double_t currentPos[3],currentSigma[3];
352 0 : vert->GetXYZ(currentPos);
353 0 : vert->GetSigmaXYZ(currentSigma);
354 : Bool_t goon = kTRUE;
355 0 : for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
356 0 : if(!goon)return;
357 0 : for(Int_t i=0;i<3;i++){
358 0 : fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
359 0 : fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
360 0 : fAverPosSum[i]+=currentPos[i];
361 : }
362 0 : for(Int_t i=0;i<3;i++){
363 0 : for(Int_t j=i;j<3;j++){
364 0 : fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
365 : }
366 : }
367 :
368 0 : fTotContributors+=vert->GetNContributors();
369 0 : fVertexXY->Fill(currentPos[0],currentPos[1]);
370 0 : fVertexZ->Fill(currentPos[2]);
371 :
372 0 : fNoEventsContr++;
373 0 : }
374 :
375 : //______________________________________________________________________
376 : Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
377 : // compute mean vertex
378 : // called once with killOutliers = kFALSE and then again with
379 : // killOutliers = kTRUE
380 :
381 : static Bool_t preliminary = kTRUE;
382 : static Int_t oldnumbevt = 0;
383 0 : if(!(preliminary || killOutliers))return kTRUE; //ComputeMean is
384 : // executed only once with argument kFALSE
385 0 : Double_t wpos[3];
386 0 : for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
387 :
388 0 : Int_t istart = oldnumbevt;
389 0 : if(killOutliers)istart = 0;
390 0 : if(killOutliers && preliminary){
391 0 : preliminary = kFALSE;
392 0 : Reset(kTRUE,kFALSE);
393 0 : }
394 0 : oldnumbevt = fVertArray.GetEntries();
395 0 : for(Int_t i=istart;i<fVertArray.GetEntries();i++){
396 : AliESDVertex* vert = NULL;
397 :
398 : // second pass
399 0 : if(killOutliers && fAccEvents.TestBitNumber(i)){
400 0 : vert=(AliESDVertex*)fVertArray[i];
401 0 : Double_t dist=(vert->GetX()-wpos[0])*(vert->GetX()-wpos[0]);
402 0 : dist+=(vert->GetY()-wpos[1])*(vert->GetY()-wpos[1]);
403 0 : dist=sqrt(dist)*10.; // distance in mm
404 0 : fDistH->Fill(dist);
405 0 : if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
406 0 : }
407 :
408 0 : if(!fAccEvents.TestBitNumber(i))continue;
409 0 : vert=(AliESDVertex*)fVertArray[i];
410 0 : AddToMean(vert);
411 0 : }
412 0 : Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
413 0 : if(bad) {
414 0 : if(killOutliers){
415 : // when the control reachs this point, the preliminary evaluation of the
416 : // diamond region has to be redone. So a general reset must be issued
417 : // if this occurs, it is likely that the cuts are badly defined
418 0 : ResetArray();
419 0 : Reset(kFALSE,kTRUE);
420 0 : preliminary = kTRUE;
421 0 : oldnumbevt = 0;
422 0 : }
423 0 : return kFALSE;
424 : }
425 0 : Double_t nevents = fNoEventsContr;
426 0 : for(Int_t i=0;i<3;i++){
427 0 : fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
428 0 : fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
429 0 : fAverPos[i] = fAverPosSum[i]/nevents;
430 : }
431 0 : for(Int_t i=0;i<3;i++){
432 0 : for(Int_t j=i;j<3;j++){
433 0 : fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
434 0 : fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
435 : }
436 : }
437 0 : if(killOutliers)ResetArray();
438 0 : fAverContributors = fTotContributors/nevents;
439 : return kTRUE;
440 0 : }
|
