Clean-up of include files
[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
21 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
22 ClassImp(AliITSMeanVertexer)
23
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 // AliITSMeanVertexer mv("RawDataFileName");
33 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
34 // ....  optional setters ....
35 // mv.Reconstruct();  // SPD local reconstruction
36 // mv.DoVertices(); 
37 // Resulting AliMeanVertex object is stored in a file named fMVFileName
38 ///////////////////////////////////////////////////////////////////////
39
40 /* $Id$ */
41
42 //______________________________________________________________________
43 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
44 fDetTypeRec(NULL),
45 fVertexXY(NULL),
46 fVertexZ(NULL),
47 fNoEventsContr(0),
48 fTotContributors(0.),
49 fAverContributors(0.), 
50 fFilterOnContributors(0),
51 fMode(mode),
52 fVertexer(NULL),
53 fAccEvents(fgkMaxNumOfEvents), 
54 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
55 fClu0(NULL),
56 fIndex(0),
57 fErrXCut(0.),
58 fRCut(0.),
59 fLowSPD0(0),
60 fHighSPD0(0)
61 {
62   // Default Constructor
63   Reset();
64
65   // Histograms initialization
66   const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
67   const Float_t xDelta = 0.02, yDelta = 0.02, 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   fClu0 = new UInt_t [fgkMaxNumOfEvents];
78   for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
79   fAccEvents.ResetAllBits(kTRUE);
80 }
81
82 //______________________________________________________________________
83 void AliITSMeanVertexer::Reset(){
84   // Reset averages 
85   for(Int_t i=0;i<3;i++){
86     fWeighPosSum[i] = 0.;
87     fWeighSigSum[i] = 0.;
88     fAverPosSum[i] = 0.;
89     fWeighPos[i] = 0.;
90     fWeighSig[i] = 0.;
91     fAverPos[i] = 0.;
92     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
93     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
94     fNoEventsContr=0;
95     fTotContributors = 0.;
96     if(fVertexXY)fVertexXY->Reset();
97     if(fVertexZ)fVertexZ->Reset();
98   }
99 }
100 //______________________________________________________________________
101 Bool_t AliITSMeanVertexer::Init() {
102   // Initialize filters
103   // Initialize geometry
104   // Initialize ITS classes
105  
106   AliGeomManager::LoadGeometry();
107   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
108
109   AliITSInitGeometry initgeom;
110   AliITSgeom *geom = initgeom.CreateAliITSgeom();
111   if (!geom) return kFALSE;
112   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
113
114   fDetTypeRec = new AliITSDetTypeRec();
115   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
116   fDetTypeRec->SetITSgeom(geom);
117   fDetTypeRec->SetDefaults();
118   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
119
120   // Initialize filter values to their defaults
121   SetFilterOnContributors();
122   SetCutOnErrX();
123   SetCutOnR();
124   SetCutOnCls();
125
126   // Instatiate vertexer
127   if (!fMode) {
128     fVertexer = new AliITSVertexer3DTapan(1000);
129   }
130   else {
131     fVertexer = new AliITSVertexer3D();
132     fVertexer->SetDetTypeRec(fDetTypeRec);
133     AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
134     alias->SetWideFiducialRegion(40.,0.5);
135     alias->SetNarrowFiducialRegion(0.5,0.5);
136     alias->SetDeltaPhiCuts(0.5,0.025);
137     alias->SetDCACut(0.1);
138     alias->SetPileupAlgo(3);
139     fVertexer->SetComputeMultiplicity(kFALSE);
140   }
141   return kTRUE;
142 }
143
144 //______________________________________________________________________
145 AliITSMeanVertexer::~AliITSMeanVertexer() {
146   // Destructor
147   delete fDetTypeRec;
148   delete fVertexXY;
149   delete fVertexZ;
150   delete fVertexer;
151 }
152
153 //______________________________________________________________________
154 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
155   // Performs SPD local reconstruction
156   // and vertex finding
157   // returns true in case a vertex is found
158
159   // Run SPD cluster finder
160   AliITSRecPointContainer::Instance()->PrepareToRead();
161   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
162   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
163
164   Bool_t vtxOK = kFALSE;
165   UInt_t nl1=0;
166   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
167   if (!fMode) {
168     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
169   }
170   else {
171     AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
172     nl1=rpcont->GetNClustersInLayerFast(1);
173     if(Filter(vtx,nl1)) vtxOK = kTRUE;
174   }
175   delete clustersTree;
176   if (vtxOK && fMode){
177     new(fVertArray[fIndex])AliESDVertex(*vtx);
178     fClu0[fIndex]=nl1;
179     fIndex++;
180   }
181   if (vtx) delete vtx;
182
183   return vtxOK;
184 }
185
186 //______________________________________________________________________
187 void AliITSMeanVertexer::WriteVertices(const char *filename){
188   // Compute mean vertex and
189   // store it along with the histograms
190   // in a file
191   
192   TFile fmv(filename,"update");
193   Bool_t itisOK = kFALSE;
194   if(ComputeMean(kFALSE)){
195     if(ComputeMean(kTRUE)){
196       Double_t cov[6];
197       cov[0] =  fAverPosSq[0][0];  // variance x
198       cov[1] =  fAverPosSq[0][1];  // cov xy
199       cov[2] =  fAverPosSq[1][1];  // variance y
200       cov[3] =  fAverPosSq[0][2];  // cov xz
201       cov[4] =  fAverPosSq[1][2];  // cov yz
202       cov[5] =  fAverPosSq[2][2];  // variance z
203       AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
204       mv.SetTitle("Mean Vertex");
205       mv.SetName("MeanVertex");
206       // we have to add chi2 here
207       AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
208
209       mv.Write(mv.GetName(),TObject::kOverwrite);
210       vtx.Write(vtx.GetName(),TObject::kOverwrite);
211       itisOK = kTRUE;
212     }
213   }
214
215   if(!itisOK){
216     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
217   }
218
219   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
220   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
221   fmv.Close();
222 }
223
224 //______________________________________________________________________
225 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
226   // Apply selection criteria to events
227   Bool_t status = kFALSE;
228   if(!vert)return status;
229   // Remove vertices reconstructed with vertexerZ
230   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
231   Int_t ncontr = vert->GetNContributors();
232   AliDebug(1,Form("Number of contributors = %d",ncontr));
233   Double_t ex = vert->GetXRes();
234   if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
235      ((ex*1000.)<fErrXCut)) {
236     status = kTRUE;
237     fAccEvents.SetBitNumber(fIndex);
238   }
239   return status;
240 }
241
242 //______________________________________________________________________
243 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
244   // update mean vertex
245   Double_t currentPos[3],currentSigma[3];
246   vert->GetXYZ(currentPos);
247   vert->GetSigmaXYZ(currentSigma);
248   Bool_t goon = kTRUE;
249   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
250   if(!goon)return;
251   for(Int_t i=0;i<3;i++){
252     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
253     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
254     fAverPosSum[i]+=currentPos[i];
255   }
256   for(Int_t i=0;i<3;i++){
257     for(Int_t j=i;j<3;j++){
258       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
259     }
260   }
261
262   fTotContributors+=vert->GetNContributors();
263   fVertexXY->Fill(currentPos[0],currentPos[1]);
264   fVertexZ->Fill(currentPos[2]);
265
266   fNoEventsContr++;
267 }
268
269 //______________________________________________________________________
270  Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
271    // compute mean vertex
272    // called once with killOutliers = kFALSE and then again with
273    // killOutliers = kTRUE
274    Double_t wpos[3];
275    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
276    Reset();
277    for(Int_t i=0;i<fVertArray.GetEntries();i++){
278      AliESDVertex* vert = NULL;
279
280      // second pass
281      if(killOutliers && fAccEvents.TestBitNumber(i)){
282        vert=(AliESDVertex*)fVertArray[i];
283        Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
284        dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
285        dist=sqrt(dist)*10.;    // distance in mm
286        if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
287      }
288      
289      if(!fAccEvents.TestBitNumber(i))continue;
290      vert=(AliESDVertex*)fVertArray[i];
291      AddToMean(vert);
292    }
293    if(fNoEventsContr < 2) return kFALSE;
294    Double_t nevents = fNoEventsContr;
295    for(Int_t i=0;i<3;i++){
296      fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; 
297      fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
298      fAverPos[i] = fAverPosSum[i]/nevents;
299    } 
300    for(Int_t i=0;i<3;i++){
301      for(Int_t j=i;j<3;j++){
302        fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
303        fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
304      }
305    } 
306    
307    fAverContributors = fTotContributors/nevents;
308    return kTRUE;
309  }