New binning for default histograms. Updated defaults (checked on PbPb run 137161)
[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 //                                                                   //
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)                             //
32 // Usage:                                                            //
33 // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm    //
34 ///////////////////////////////////////////////////////////////////////
35
36 /* $Id$ */
37
38 //______________________________________________________________________
39 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
40 fDetTypeRec(NULL),
41 fVertexXY(NULL),
42 fVertexZ(NULL),
43 fNoEventsContr(0),
44 fTotContributors(0.),
45 fAverContributors(0.), 
46 fFilterOnContributors(0),
47 fMode(mode),
48 fVertexer(NULL),
49 fAccEvents(fgkMaxNumOfEvents), 
50 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
51 fClu0(NULL),
52 fIndex(0),
53 fErrXCut(0.),
54 fRCut(0.),
55 fLowSPD0(0),
56 fHighSPD0(0),
57 fMultH(NULL),
58 fErrXH(NULL),
59 fMultHa(NULL),
60 fErrXHa(NULL),
61 fDistH(NULL)
62 {
63   // Default Constructor
64   Reset();
65
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);
88 }
89
90 //______________________________________________________________________
91 void AliITSMeanVertexer::Reset(){
92   // Reset averages 
93   for(Int_t i=0;i<3;i++){
94     fWeighPosSum[i] = 0.;
95     fWeighSigSum[i] = 0.;
96     fAverPosSum[i] = 0.;
97     fWeighPos[i] = 0.;
98     fWeighSig[i] = 0.;
99     fAverPos[i] = 0.;
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.;
102     fNoEventsContr=0;
103     fTotContributors = 0.;
104     if(fVertexXY)fVertexXY->Reset();
105     if(fVertexZ)fVertexZ->Reset();
106   }
107 }
108 //______________________________________________________________________
109 Bool_t AliITSMeanVertexer::Init() {
110   // Initialize filters
111   // Initialize geometry
112   // Initialize ITS classes
113  
114   AliGeomManager::LoadGeometry();
115   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
116
117   AliITSInitGeometry initgeom;
118   AliITSgeom *geom = initgeom.CreateAliITSgeom();
119   if (!geom) return kFALSE;
120   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
121
122   fDetTypeRec = new AliITSDetTypeRec();
123   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
124   fDetTypeRec->SetITSgeom(geom);
125   fDetTypeRec->SetDefaults();
126   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
127
128   // Initialize filter values to their defaults
129   SetFilterOnContributors();
130   SetCutOnErrX();
131   SetCutOnR();
132   SetCutOnCls();
133
134   // Instatiate vertexer
135   if (!fMode) {
136     fVertexer = new AliITSVertexer3DTapan(1000);
137   }
138   else {
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);
148   }
149   return kTRUE;
150 }
151
152 //______________________________________________________________________
153 AliITSMeanVertexer::~AliITSMeanVertexer() {
154   // Destructor
155   delete fDetTypeRec;
156   delete fVertexXY;
157   delete fVertexZ;
158   delete fMultH;
159   delete fErrXH;
160   delete fMultHa;
161   delete fErrXHa;
162   delete fDistH;
163   delete fVertexer;
164 }
165
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
171
172   // Run SPD cluster finder
173   static Int_t evcount = -1;
174   if(evcount <0){
175     evcount++;
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));
178   }
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");
183
184   Bool_t vtxOK = kFALSE;
185   UInt_t nl1=0;
186   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
187   if (!fMode) {
188     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
189   }
190   else {
191     AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
192     nl1=rpcont->GetNClustersInLayerFast(1);
193     if(Filter(vtx,nl1)) vtxOK = kTRUE;
194     /*    if(vtx){
195       if(vtxOK){
196         printf("The vertex is OK\n");
197       }
198       else {
199         printf("The vertex is NOT OK\n");
200       }
201       vtx->PrintStatus();
202     }
203     else {
204       printf("The vertex was not reconstructed\n");
205       }  */
206   }
207   delete clustersTree;
208   if (vtxOK && fMode){
209     new(fVertArray[fIndex])AliESDVertex(*vtx);
210     fClu0[fIndex]=nl1;
211     fAccEvents.SetBitNumber(fIndex);
212     fIndex++;
213     if(fIndex>=fgkMaxNumOfEvents){
214       if(ComputeMean(kFALSE)){
215         if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
216       }
217     }
218   }
219   if (vtx) delete vtx;
220
221   return vtxOK;
222 }
223
224 //______________________________________________________________________
225 void AliITSMeanVertexer::WriteVertices(const char *filename){
226   // Compute mean vertex and
227   // store it along with the histograms
228   // in a file
229   
230   TFile fmv(filename,"update");
231   Bool_t itisOK = kFALSE;
232   if(ComputeMean(kFALSE)){
233     if(ComputeMean(kTRUE)){
234       Double_t cov[6];
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");
246
247       mv.Write(mv.GetName(),TObject::kOverwrite);
248       vtx.Write(vtx.GetName(),TObject::kOverwrite);
249       itisOK = kTRUE;
250     }
251   }
252
253   if(!itisOK){
254     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
255   }
256
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);
264   fmv.Close();
265 }
266
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();
277   fMultH->Fill(mult);
278   fErrXH->Fill(ex*1000.);
279   if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
280      ((ex*1000.)<fErrXCut)) {
281     status = kTRUE;
282     fMultHa->Fill(mult);
283     fErrXHa->Fill(ex*1000.);
284   }
285   return status;
286 }
287
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);
294   Bool_t goon = kTRUE;
295   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
296   if(!goon)return;
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];
301   }
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];
305     }
306   }
307
308   fTotContributors+=vert->GetNContributors();
309   fVertexXY->Fill(currentPos[0],currentPos[1]);
310   fVertexZ->Fill(currentPos[2]);
311
312   fNoEventsContr++;
313 }
314
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
320
321    static Bool_t firstime = kTRUE;
322    if(fNoEventsContr>2 && !killOutliers)return kTRUE;  //ComputeMean is
323                              // executed only once with argument kFALSE
324    Double_t wpos[3];
325    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
326
327    if(killOutliers && firstime){
328      firstime = kFALSE;
329      Reset();
330    }
331
332    for(Int_t i=0;i<fVertArray.GetEntries();i++){
333      AliESDVertex* vert = NULL;
334
335      // second pass
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
341        fDistH->Fill(dist);
342        if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
343      }
344      
345      if(!fAccEvents.TestBitNumber(i))continue;
346      vert=(AliESDVertex*)fVertArray[i];
347      AddToMean(vert);
348    }
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;
355    } 
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];
360      }
361    } 
362    if(killOutliers)ResetArray();
363    fAverContributors = fTotContributors/nevents;
364    return kTRUE;
365  }