Changes to use the mean vertexer DA
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
1 #include <TFile.h>
2 #include <TH1.h>
3 #include <TH2.h>
4 #include "AliGeomManager.h"
5 #include "AliITSDetTypeRec.h"
6 #include "AliITSInitGeometry.h"
7 #include "AliITSMeanVertexer.h"
8 #include "AliITSRecPointContainer.h"
9 #include "AliITSLoader.h"
10 #include "AliLog.h"
11 #include "AliRawReader.h"
12 #include "AliRawReaderDate.h"
13 #include "AliRawReaderRoot.h"
14 #include "AliRunLoader.h"
15 #include "AliITSVertexer3D.h"
16 #include "AliITSVertexer3DTapan.h"
17 #include "AliESDVertex.h"
18 #include "AliMeanVertex.h"
19 #include "AliMultiplicity.h"
20
21 ClassImp(AliITSMeanVertexer)
22
23 ///////////////////////////////////////////////////////////////////////
24 //                                                                   //
25 // Class to compute vertex position using SPD local reconstruction   //
26 // An average vertex position using all the events                   //
27 // is built and saved                                                //
28 // Individual vertices can be optionally stored                      //
29 // Origin: M.Masera  (masera@to.infn.it)                             //
30 // Usage:
31 // AliITSMeanVertexer mv("RawDataFileName");
32 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
33 // ....  optional setters ....
34 // mv.Reconstruct();  // SPD local reconstruction
35 // mv.DoVertices(); 
36 // Resulting AliMeanVertex object is stored in a file named fMVFileName
37 ///////////////////////////////////////////////////////////////////////
38
39 /* $Id$ */
40
41 //______________________________________________________________________
42 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
43 fDetTypeRec(NULL),
44 fVertexXY(NULL),
45 fVertexZ(NULL),
46 fNoEventsContr(0),
47 fTotTracklets(0.),
48 fAverTracklets(0.),
49 fTotTrackletsSq(0.),
50 fSigmaOnAverTracks(0.), 
51 fFilterOnContributors(0),
52 fFilterOnTracklets(0),
53 fMode(mode),
54 fVertexer(NULL)
55 {
56   // Default Constructor
57   for(Int_t i=0;i<3;i++){
58     fWeighPosSum[i] = 0.;
59     fWeighSigSum[i] = 0.;
60     fAverPosSum[i] = 0.;
61     fWeighPos[i] = 0.;
62     fWeighSig[i] = 0.;
63     fAverPos[i] = 0.;
64     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
65     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
66   }
67
68   // Histograms initialization
69   const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
70   const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
71   fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
72                        2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
73                        2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
74   fVertexXY->SetXTitle("X , cm");
75   fVertexXY->SetYTitle("Y , cm");
76   fVertexXY->SetOption("colz");
77   fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
78                        2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
79   fVertexZ->SetXTitle("Z , cm");
80 }
81
82 //______________________________________________________________________
83 Bool_t AliITSMeanVertexer::Init() {
84   // Initialize filters
85   // Initialize geometry
86   // Initialize ITS classes
87  
88   AliGeomManager::LoadGeometry();
89   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
90
91   AliITSInitGeometry initgeom;
92   AliITSgeom *geom = initgeom.CreateAliITSgeom();
93   if (!geom) return kFALSE;
94   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
95
96   fDetTypeRec = new AliITSDetTypeRec();
97   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
98   fDetTypeRec->SetITSgeom(geom);
99   fDetTypeRec->SetDefaults();
100   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
101
102   // Initialize filter values to their defaults
103   SetFilterOnContributors();
104   SetFilterOnTracklets();
105
106   // Instatiate vertexer
107   if (!fMode) {
108     fVertexer = new AliITSVertexer3DTapan(1000);
109   }
110   else {
111     fVertexer = new AliITSVertexer3D();
112     fVertexer->SetDetTypeRec(fDetTypeRec);
113     fVertexer->SetComputeMultiplicity(kTRUE);
114   }
115   return kTRUE;
116 }
117
118 //______________________________________________________________________
119 AliITSMeanVertexer::~AliITSMeanVertexer() {
120   // Destructor
121   delete fDetTypeRec;
122   delete fVertexXY;
123   delete fVertexZ;
124   delete fVertexer;
125 }
126
127 //______________________________________________________________________
128 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
129   // Performs SPD local reconstruction
130   // and vertex finding
131   // returns true in case a vertex is found
132
133   // Run SPD cluster finder
134   AliITSRecPointContainer::Instance()->PrepareToRead();
135   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
136   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
137
138   Bool_t vtxOK = kFALSE;
139   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
140   if (!fMode) {
141     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
142   }
143   else {
144     AliMultiplicity *mult = fVertexer->GetMultiplicity();
145     if(Filter(vtx,mult)) vtxOK = kTRUE;
146   }
147   delete clustersTree;
148   if (vtxOK) AddToMean(vtx);
149   if (vtx) delete vtx;
150
151   return vtxOK;
152 }
153
154 //______________________________________________________________________
155 void AliITSMeanVertexer::WriteVertices(const char *filename){
156   // Compute mean vertex and
157   // store it along with the histograms
158   // in a file
159   
160   TFile fmv(filename,"update");
161
162   if(ComputeMean()){
163     Double_t cov[6];
164     cov[0] =  fAverPosSq[0][0];  // variance x
165     cov[1] =  fAverPosSq[0][1];  // cov xy
166     cov[2] =  fAverPosSq[1][1];  // variance y
167     cov[3] =  fAverPosSq[0][2];  // cov xz
168     cov[4] =  fAverPosSq[1][2];  // cov yz
169     cov[5] =  fAverPosSq[2][2];  // variance z
170     AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
171     mv.SetTitle("Mean Vertex");
172     mv.SetName("MeanVertex");
173     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
174     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
175     // we have to add chi2 here
176     AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
177
178     mv.Write(mv.GetName(),TObject::kOverwrite);
179     vtx.Write(vtx.GetName(),TObject::kOverwrite);
180   }
181   else {
182     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
183   }
184
185   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
186   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
187   fmv.Close();
188 }
189
190 //______________________________________________________________________
191 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
192   // Apply selection criteria to events
193   Bool_t status = kFALSE;
194   if(!vert || !mult)return status;
195   // Remove vertices reconstructed with vertexerZ
196   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
197   Int_t ncontr = vert->GetNContributors();
198   Int_t ntracklets = mult->GetNumberOfTracklets();
199   AliDebug(1,Form("Number of contributors = %d",ncontr));
200   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
201   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
202   fTotTracklets += ntracklets;
203   fTotTrackletsSq += ntracklets*ntracklets;
204   return status;
205 }
206
207 //______________________________________________________________________
208 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
209   // update mean vertex
210   Double_t currentPos[3],currentSigma[3];
211   vert->GetXYZ(currentPos);
212   vert->GetSigmaXYZ(currentSigma);
213   Bool_t goon = kTRUE;
214   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
215   if(!goon)return;
216   for(Int_t i=0;i<3;i++){
217     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
218     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
219     fAverPosSum[i]+=currentPos[i];
220   }
221   for(Int_t i=0;i<3;i++){
222     for(Int_t j=i;j<3;j++){
223       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
224     }
225   }
226
227   fVertexXY->Fill(currentPos[0],currentPos[1]);
228   fVertexZ->Fill(currentPos[2]);
229
230   fNoEventsContr++;
231 }
232
233 //______________________________________________________________________
234 Bool_t AliITSMeanVertexer::ComputeMean(){
235   // compute mean vertex
236   if(fNoEventsContr < 2) return kFALSE;
237   Double_t nevents = fNoEventsContr;
238   for(Int_t i=0;i<3;i++){
239     fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; 
240     fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
241     fAverPos[i] = fAverPosSum[i]/nevents;
242   } 
243   for(Int_t i=0;i<3;i++){
244     for(Int_t j=i;j<3;j++){
245       fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
246       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
247     }
248   } 
249
250   fAverTracklets = fTotTracklets/nevents;
251   fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
252   fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
253   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
254   return kTRUE;
255 }