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