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