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