store also difference in local Y
[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 fZCutDiamond(0.),
53 fLowSPD0(0),
54 fHighSPD0(0),
55 fMultH(NULL),
56 fErrXH(NULL),
57 fMultHa(NULL),
58 fErrXHa(NULL),
59 fDistH(NULL),
60 fContrH(NULL),
61 fContrHa(NULL)
62 {
63   // Default Constructor
64   SetZFiducialRegion();   //  sets fZCutDiamond to its default value
65   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   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);
89
90   fClu0 = new UInt_t [fgkMaxNumOfEvents];
91   for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
92   fAccEvents.ResetAllBits(kTRUE);
93 }
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   if(fVertexXY && redefine2D){
106     counter++;
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();
111     delete fVertexXY;
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);
121   }
122   else if(fVertexXY && !redefine2D){
123     fVertexXY->Reset();
124   }
125   if(fVertexZ){
126     fVertexZ->Reset();
127     fDistH->Reset();
128     if(complete){
129       fErrXH->Reset();
130       fMultH->Reset();
131       fErrXHa->Reset(); 
132       fMultHa->Reset();
133       fContrH->Reset();
134       fContrHa->Reset();
135     }
136   }
137   for(Int_t i=0;i<3;i++){
138     fWeighPosSum[i] = 0.;
139     fWeighSigSum[i] = 0.;
140     fAverPosSum[i] = 0.;
141     fWeighPos[i] = 0.;
142     fWeighSig[i] = 0.;
143     fAverPos[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.;
146     fNoEventsContr=0;
147     fTotContributors = 0.;
148   }
149 }
150 //______________________________________________________________________
151 Bool_t AliITSMeanVertexer::Init() {
152   // Initialize filters
153   // Initialize geometry
154   // Initialize ITS classes
155  
156   AliGeomManager::LoadGeometry();
157   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
158
159   AliITSInitGeometry initgeom;
160   AliITSgeom *geom = initgeom.CreateAliITSgeom();
161   if (!geom) return kFALSE;
162   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
163
164   fDetTypeRec = new AliITSDetTypeRec();
165   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
166   fDetTypeRec->SetITSgeom(geom);
167   fDetTypeRec->SetDefaults();
168   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
169
170   // Initialize filter values to their defaults
171   SetFilterOnContributors();
172   SetCutOnErrX();
173   SetCutOnR();
174   SetCutOnCls();
175
176   // Instatiate vertexer
177   if (!fMode) {
178     fVertexer = new AliITSVertexer3DTapan(1000);
179   }
180   else {
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);
190   }
191   return kTRUE;
192 }
193
194 //______________________________________________________________________
195 AliITSMeanVertexer::~AliITSMeanVertexer() {
196   // Destructor
197   delete fDetTypeRec;
198   delete fVertexXY;
199   delete fVertexZ;
200   delete fMultH;
201   delete fErrXH;
202   delete fMultHa;
203   delete fErrXHa;
204   delete fDistH;
205   delete fContrH;
206   delete fContrHa;
207   delete fVertexer;
208 }
209
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
215
216   // Run SPD cluster finder
217   static Int_t evcount = -1;
218   if(evcount <0){
219     evcount++;
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));
222   }
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");
227
228   Bool_t vtxOK = kFALSE;
229   UInt_t nl1=0;
230   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
231   if (!fMode) {
232     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
233   }
234   else {
235     AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
236     nl1=rpcont->GetNClustersInLayerFast(1);
237     if(Filter(vtx,nl1)) vtxOK = kTRUE;
238     /*    if(vtx){
239       if(vtxOK){
240         printf("The vertex is OK\n");
241       }
242       else {
243         printf("The vertex is NOT OK\n");
244       }
245       vtx->PrintStatus();
246     }
247     else {
248       printf("The vertex was not reconstructed\n");
249       }  */
250   }
251   delete clustersTree;
252   if (vtxOK && fMode){
253     new(fVertArray[fIndex])AliESDVertex(*vtx);
254     fClu0[fIndex]=nl1;
255     fAccEvents.SetBitNumber(fIndex);
256     fIndex++;
257     if(fIndex>=fgkMaxNumOfEvents){
258       if(ComputeMean(kFALSE)){
259         if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
260       }
261     }
262   }
263   if (vtx) delete vtx;
264
265   return vtxOK;
266 }
267
268 //______________________________________________________________________
269 void AliITSMeanVertexer::WriteVertices(const char *filename){
270   // Compute mean vertex and
271   // store it along with the histograms
272   // in a file
273   
274   TFile fmv(filename,"update");
275   Bool_t itisOK = kFALSE;
276   if(ComputeMean(kFALSE)){
277     if(ComputeMean(kTRUE)){
278       Double_t cov[6];
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");
295
296       mv.Write(mv.GetName(),TObject::kOverwrite);
297       vtx.Write(vtx.GetName(),TObject::kOverwrite);
298       itisOK = kTRUE;
299     }
300   }
301
302   if(!itisOK){
303     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
304   }
305
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);
315   fmv.Close();
316 }
317
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();
328   fMultH->Fill(mult);
329   fErrXH->Fill(ex*1000.);
330   fContrH->Fill((Float_t)ncontr);
331   if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
332      ((ex*1000.)<fErrXCut)) {
333     status = kTRUE;
334     fMultHa->Fill(mult);
335     fErrXHa->Fill(ex*1000.);
336     fContrHa->Fill((Float_t)ncontr);
337   }
338   return status;
339 }
340
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);
347   Bool_t goon = kTRUE;
348   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
349   if(!goon)return;
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];
354   }
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];
358     }
359   }
360
361   fTotContributors+=vert->GetNContributors();
362   fVertexXY->Fill(currentPos[0],currentPos[1]);
363   fVertexZ->Fill(currentPos[2]);
364
365   fNoEventsContr++;
366 }
367
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
373
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
378    Double_t wpos[3];
379    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
380
381    Int_t istart = oldnumbevt;
382    if(killOutliers)istart = 0;
383    if(killOutliers && preliminary){
384      preliminary = kFALSE;
385      Reset(kTRUE,kFALSE);
386    }
387    oldnumbevt = fVertArray.GetEntries();
388    for(Int_t i=istart;i<fVertArray.GetEntries();i++){
389      AliESDVertex* vert = NULL;
390
391      // second pass
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
397        fDistH->Fill(dist);
398        if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
399      }
400      
401      if(!fAccEvents.TestBitNumber(i))continue;
402      vert=(AliESDVertex*)fVertArray[i];
403      AddToMean(vert);
404    }
405    Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
406    if(bad) {
407      if(killOutliers){  
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
411        ResetArray();
412        Reset(kFALSE,kTRUE);
413        preliminary = kTRUE;
414        oldnumbevt = 0;
415      }
416      return kFALSE;
417    }
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;
423    } 
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];
428      }
429    } 
430    if(killOutliers)ResetArray();
431    fAverContributors = fTotContributors/nevents;
432    return kTRUE;
433  }