]>
Commit | Line | Data |
---|---|---|
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 |
8d7f6205 | 66 | const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0; |
37add1d1 | 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); | |
63845149 | 71 | fVertexXY->SetXTitle("X , cm"); |
72 | fVertexXY->SetYTitle("Y , cm"); | |
73 | fVertexXY->SetOption("colz"); | |
37add1d1 | 74 | fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile", |
75 | 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit); | |
63845149 | 76 | fVertexZ->SetXTitle("Z , cm"); |
a3a3a28f | 77 | } |
78 | ||
79 | //______________________________________________________________________ | |
37add1d1 | 80 | Bool_t AliITSMeanVertexer::Init() { |
81 | // Initialize filters | |
a3a3a28f | 82 | // Initialize geometry |
37add1d1 | 83 | // Initialize ITS classes |
a3a3a28f | 84 | |
37add1d1 | 85 | AliGeomManager::LoadGeometry(); |
86 | if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE; | |
a3a3a28f | 87 | |
88 | AliITSInitGeometry initgeom; | |
89 | AliITSgeom *geom = initgeom.CreateAliITSgeom(); | |
37add1d1 | 90 | if (!geom) return kFALSE; |
a3a3a28f | 91 | printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data()); |
37add1d1 | 92 | |
cf0b1cea | 93 | fDetTypeRec = new AliITSDetTypeRec(); |
762133f4 | 94 | fDetTypeRec->SetLoadOnlySPDCalib(kTRUE); |
37add1d1 | 95 | fDetTypeRec->SetITSgeom(geom); |
96 | fDetTypeRec->SetDefaults(); | |
97 | fDetTypeRec->SetDefaultClusterFindersV2(kTRUE); | |
98 | ||
a3a3a28f | 99 | // Initialize filter values to their defaults |
100 | SetFilterOnContributors(); | |
101 | SetFilterOnTracklets(); | |
a3a3a28f | 102 | |
37add1d1 | 103 | return kTRUE; |
a3a3a28f | 104 | } |
105 | ||
106 | //______________________________________________________________________ | |
107 | AliITSMeanVertexer::~AliITSMeanVertexer() { | |
108 | // Destructor | |
37add1d1 | 109 | delete fDetTypeRec; |
110 | delete fVertexXY; | |
111 | delete fVertexZ; | |
a3a3a28f | 112 | } |
113 | ||
114 | //______________________________________________________________________ | |
37add1d1 | 115 | Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){ |
a3a3a28f | 116 | // Performs SPD local reconstruction |
37add1d1 | 117 | // and vertex finding |
118 | // returns true in case a vertex is found | |
119 | ||
120 | // Run SPD cluster finder | |
121 | TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree | |
122 | fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD"); | |
123 | ||
124 | Bool_t vtxOK = kFALSE; | |
125 | AliESDVertex *vtx = NULL; | |
126 | // Run Tapan's vertex-finder | |
127 | if (!mode) { | |
128 | AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000); | |
129 | vtx = vertexer1->FindVertexForCurrentEvent(clustersTree); | |
130 | delete vertexer1; | |
131 | if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE; | |
a3a3a28f | 132 | } |
37add1d1 | 133 | else { |
134 | // Run standard vertexer3d | |
135 | AliITSVertexer3D *vertexer2 = new AliITSVertexer3D(); | |
307ebfc0 | 136 | vertexer2->SetDetTypeRec(fDetTypeRec); |
cf0b1cea | 137 | vtx = vertexer2->FindVertexForCurrentEvent(clustersTree); |
37add1d1 | 138 | AliMultiplicity *mult = vertexer2->GetMultiplicity(); |
139 | delete vertexer2; | |
140 | if(Filter(vtx,mult)) vtxOK = kTRUE; | |
a3a3a28f | 141 | } |
37add1d1 | 142 | delete clustersTree; |
143 | if (vtxOK) AddToMean(vtx); | |
144 | if (vtx) delete vtx; | |
a3a3a28f | 145 | |
37add1d1 | 146 | return vtxOK; |
a3a3a28f | 147 | } |
148 | ||
149 | //______________________________________________________________________ | |
37add1d1 | 150 | void AliITSMeanVertexer::WriteVertices(const char *filename){ |
151 | // Compute mean vertex and | |
152 | // store it along with the histograms | |
153 | // in a file | |
154 | ||
155 | TFile fmv(filename,"update"); | |
156 | ||
a3a3a28f | 157 | if(ComputeMean()){ |
158 | Double_t cov[6]; | |
159 | cov[0] = fAverPosSq[0][0]; // variance x | |
160 | cov[1] = fAverPosSq[0][1]; // cov xy | |
161 | cov[2] = fAverPosSq[1][1]; // variance y | |
162 | cov[3] = fAverPosSq[0][2]; // cov xz | |
163 | cov[4] = fAverPosSq[1][2]; // cov yz | |
164 | cov[5] = fAverPosSq[2][2]; // variance z | |
37add1d1 | 165 | AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks); |
a3a3a28f | 166 | mv.SetTitle("Mean Vertex"); |
37add1d1 | 167 | mv.SetName("MeanVertex"); |
a3a3a28f | 168 | AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets())); |
169 | AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks())); | |
37add1d1 | 170 | // we have to add chi2 here |
762133f4 | 171 | AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos"); |
37add1d1 | 172 | |
173 | mv.Write(mv.GetName(),TObject::kOverwrite); | |
174 | vtx.Write(vtx.GetName(),TObject::kOverwrite); | |
a3a3a28f | 175 | } |
176 | else { | |
177 | AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr)); | |
178 | } | |
37add1d1 | 179 | |
180 | fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite); | |
181 | fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite); | |
182 | fmv.Close(); | |
a3a3a28f | 183 | } |
184 | ||
185 | //______________________________________________________________________ | |
186 | Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){ | |
187 | // Apply selection criteria to events | |
188 | Bool_t status = kFALSE; | |
189 | if(!vert || !mult)return status; | |
8d7f6205 | 190 | // Remove vertices reconstructed with vertexerZ |
191 | if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status; | |
a3a3a28f | 192 | Int_t ncontr = vert->GetNContributors(); |
193 | Int_t ntracklets = mult->GetNumberOfTracklets(); | |
194 | AliDebug(1,Form("Number of contributors = %d",ncontr)); | |
195 | AliDebug(1,Form("Number of tracklets = %d",ntracklets)); | |
196 | if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE; | |
762133f4 | 197 | fTotTracklets += ntracklets; |
198 | fTotTrackletsSq += ntracklets*ntracklets; | |
a3a3a28f | 199 | return status; |
200 | } | |
201 | ||
202 | //______________________________________________________________________ | |
203 | void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){ | |
204 | // update mean vertex | |
205 | Double_t currentPos[3],currentSigma[3]; | |
206 | vert->GetXYZ(currentPos); | |
207 | vert->GetSigmaXYZ(currentSigma); | |
208 | Bool_t goon = kTRUE; | |
209 | for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE; | |
210 | if(!goon)return; | |
211 | for(Int_t i=0;i<3;i++){ | |
762133f4 | 212 | fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i]; |
213 | fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i]; | |
214 | fAverPosSum[i]+=currentPos[i]; | |
a3a3a28f | 215 | } |
216 | for(Int_t i=0;i<3;i++){ | |
217 | for(Int_t j=i;j<3;j++){ | |
762133f4 | 218 | fAverPosSqSum[i][j] += currentPos[i] * currentPos[j]; |
a3a3a28f | 219 | } |
220 | } | |
37add1d1 | 221 | |
222 | fVertexXY->Fill(currentPos[0],currentPos[1]); | |
223 | fVertexZ->Fill(currentPos[2]); | |
224 | ||
a3a3a28f | 225 | fNoEventsContr++; |
226 | } | |
227 | ||
228 | //______________________________________________________________________ | |
229 | Bool_t AliITSMeanVertexer::ComputeMean(){ | |
230 | // compute mean vertex | |
231 | if(fNoEventsContr < 2) return kFALSE; | |
232 | Double_t nevents = fNoEventsContr; | |
233 | for(Int_t i=0;i<3;i++){ | |
762133f4 | 234 | fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; |
235 | fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]); | |
236 | fAverPos[i] = fAverPosSum[i]/nevents; | |
a3a3a28f | 237 | } |
238 | for(Int_t i=0;i<3;i++){ | |
239 | for(Int_t j=i;j<3;j++){ | |
762133f4 | 240 | fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.); |
a3a3a28f | 241 | fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j]; |
242 | } | |
243 | } | |
762133f4 | 244 | |
245 | fAverTracklets = fTotTracklets/nevents; | |
246 | fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1); | |
247 | fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets; | |
a3a3a28f | 248 | fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks); |
249 | return kTRUE; | |
250 | } |