4 #include "AliGeomManager.h"
5 #include "AliITSDetTypeRec.h"
6 #include "AliITSInitGeometry.h"
7 #include "AliITSMeanVertexer.h"
8 #include "AliITSLoader.h"
10 #include "AliRawReader.h"
11 #include "AliRawReaderDate.h"
12 #include "AliRawReaderRoot.h"
13 #include "AliRunLoader.h"
14 #include "AliITSVertexer3D.h"
15 #include "AliITSVertexer3DTapan.h"
16 #include "AliESDVertex.h"
17 #include "AliMeanVertex.h"
18 #include "AliMultiplicity.h"
20 ClassImp(AliITSMeanVertexer)
22 ///////////////////////////////////////////////////////////////////////
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) //
30 // AliITSMeanVertexer mv("RawDataFileName");
31 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
32 // .... optional setters ....
33 // mv.Reconstruct(); // SPD local reconstruction
35 // Resulting AliMeanVertex object is stored in a file named fMVFileName
36 ///////////////////////////////////////////////////////////////////////
40 //______________________________________________________________________
41 AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
49 fSigmaOnAverTracks(0.),
50 fFilterOnContributors(0),
53 // Default Constructor
54 for(Int_t i=0;i<3;i++){
61 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
62 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
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);
75 //______________________________________________________________________
76 Bool_t AliITSMeanVertexer::Init() {
78 // Initialize geometry
79 // Initialize ITS classes
81 AliGeomManager::LoadGeometry();
82 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
84 AliITSInitGeometry initgeom;
85 AliITSgeom *geom = initgeom.CreateAliITSgeom();
86 if (!geom) return kFALSE;
87 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
89 fDetTypeRec = new AliITSDetTypeRec();
90 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
91 fDetTypeRec->SetITSgeom(geom);
92 fDetTypeRec->SetDefaults();
93 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
95 // Initialize filter values to their defaults
96 SetFilterOnContributors();
97 SetFilterOnTracklets();
102 //______________________________________________________________________
103 AliITSMeanVertexer::~AliITSMeanVertexer() {
110 //______________________________________________________________________
111 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
112 // Performs SPD local reconstruction
113 // and vertex finding
114 // returns true in case a vertex is found
116 // Run SPD cluster finder
117 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
118 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
120 Bool_t vtxOK = kFALSE;
121 AliESDVertex *vtx = NULL;
122 // Run Tapan's vertex-finder
124 AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
125 vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
127 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
130 // Run standard vertexer3d
131 AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
132 vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
133 AliMultiplicity *mult = vertexer2->GetMultiplicity();
135 if(Filter(vtx,mult)) vtxOK = kTRUE;
138 if (vtxOK) AddToMean(vtx);
144 //______________________________________________________________________
145 void AliITSMeanVertexer::WriteVertices(const char *filename){
146 // Compute mean vertex and
147 // store it along with the histograms
150 TFile fmv(filename,"update");
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
160 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
161 mv.SetTitle("Mean Vertex");
162 mv.SetName("MeanVertex");
163 AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
164 AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
165 // we have to add chi2 here
166 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
168 mv.Write(mv.GetName(),TObject::kOverwrite);
169 vtx.Write(vtx.GetName(),TObject::kOverwrite);
172 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
175 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
176 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
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;
190 fTotTracklets += ntracklets;
191 fTotTrackletsSq += ntracklets*ntracklets;
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);
202 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
204 for(Int_t i=0;i<3;i++){
205 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
206 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
207 fAverPosSum[i]+=currentPos[i];
209 for(Int_t i=0;i<3;i++){
210 for(Int_t j=i;j<3;j++){
211 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
215 fVertexXY->Fill(currentPos[0],currentPos[1]);
216 fVertexZ->Fill(currentPos[2]);
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++){
227 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
228 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
229 fAverPos[i] = fAverPosSum[i]/nevents;
231 for(Int_t i=0;i<3;i++){
232 for(Int_t j=i;j<3;j++){
233 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
234 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
238 fAverTracklets = fTotTracklets/nevents;
239 fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
240 fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
241 fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);