4 #include "AliGeomManager.h"
5 #include "AliITSDetTypeRec.h"
6 #include "AliITSInitGeometry.h"
7 #include "AliITSMeanVertexer.h"
8 #include "AliITSRecPointContainer.h"
9 #include "AliITSLoader.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 #include "AliMultiplicity.h"
21 ClassImp(AliITSMeanVertexer)
23 ///////////////////////////////////////////////////////////////////////
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) //
31 // AliITSMeanVertexer mv("RawDataFileName");
32 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
33 // .... optional setters ....
34 // mv.Reconstruct(); // SPD local reconstruction
36 // Resulting AliMeanVertex object is stored in a file named fMVFileName
37 ///////////////////////////////////////////////////////////////////////
41 //______________________________________________________________________
42 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
50 fSigmaOnAverTracks(0.),
51 fFilterOnContributors(0),
52 fFilterOnTracklets(0),
56 // Default Constructor
57 for(Int_t i=0;i<3;i++){
64 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
65 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
68 // Histograms initialization
69 const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
70 const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
71 fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
72 2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
73 2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
74 fVertexXY->SetXTitle("X , cm");
75 fVertexXY->SetYTitle("Y , cm");
76 fVertexXY->SetOption("colz");
77 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
78 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
79 fVertexZ->SetXTitle("Z , cm");
82 //______________________________________________________________________
83 Bool_t AliITSMeanVertexer::Init() {
85 // Initialize geometry
86 // Initialize ITS classes
88 AliGeomManager::LoadGeometry();
89 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
91 AliITSInitGeometry initgeom;
92 AliITSgeom *geom = initgeom.CreateAliITSgeom();
93 if (!geom) return kFALSE;
94 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
96 fDetTypeRec = new AliITSDetTypeRec();
97 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
98 fDetTypeRec->SetITSgeom(geom);
99 fDetTypeRec->SetDefaults();
100 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
102 // Initialize filter values to their defaults
103 SetFilterOnContributors();
104 SetFilterOnTracklets();
106 // Instatiate vertexer
108 fVertexer = new AliITSVertexer3DTapan(1000);
111 fVertexer = new AliITSVertexer3D();
112 fVertexer->SetDetTypeRec(fDetTypeRec);
113 fVertexer->SetComputeMultiplicity(kTRUE);
118 //______________________________________________________________________
119 AliITSMeanVertexer::~AliITSMeanVertexer() {
127 //______________________________________________________________________
128 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
129 // Performs SPD local reconstruction
130 // and vertex finding
131 // returns true in case a vertex is found
133 // Run SPD cluster finder
134 AliITSRecPointContainer::Instance()->PrepareToRead();
135 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
136 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
138 Bool_t vtxOK = kFALSE;
139 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
141 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
144 AliMultiplicity *mult = fVertexer->GetMultiplicity();
145 if(Filter(vtx,mult)) vtxOK = kTRUE;
148 if (vtxOK) AddToMean(vtx);
154 //______________________________________________________________________
155 void AliITSMeanVertexer::WriteVertices(const char *filename){
156 // Compute mean vertex and
157 // store it along with the histograms
160 TFile fmv(filename,"update");
164 cov[0] = fAverPosSq[0][0]; // variance x
165 cov[1] = fAverPosSq[0][1]; // cov xy
166 cov[2] = fAverPosSq[1][1]; // variance y
167 cov[3] = fAverPosSq[0][2]; // cov xz
168 cov[4] = fAverPosSq[1][2]; // cov yz
169 cov[5] = fAverPosSq[2][2]; // variance z
170 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
171 mv.SetTitle("Mean Vertex");
172 mv.SetName("MeanVertex");
173 AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
174 AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
175 // we have to add chi2 here
176 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
178 mv.Write(mv.GetName(),TObject::kOverwrite);
179 vtx.Write(vtx.GetName(),TObject::kOverwrite);
182 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
185 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
186 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
190 //______________________________________________________________________
191 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
192 // Apply selection criteria to events
193 Bool_t status = kFALSE;
194 if(!vert || !mult)return status;
195 // Remove vertices reconstructed with vertexerZ
196 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
197 Int_t ncontr = vert->GetNContributors();
198 Int_t ntracklets = mult->GetNumberOfTracklets();
199 AliDebug(1,Form("Number of contributors = %d",ncontr));
200 AliDebug(1,Form("Number of tracklets = %d",ntracklets));
201 if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
202 fTotTracklets += ntracklets;
203 fTotTrackletsSq += ntracklets*ntracklets;
207 //______________________________________________________________________
208 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
209 // update mean vertex
210 Double_t currentPos[3],currentSigma[3];
211 vert->GetXYZ(currentPos);
212 vert->GetSigmaXYZ(currentSigma);
214 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
216 for(Int_t i=0;i<3;i++){
217 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
218 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
219 fAverPosSum[i]+=currentPos[i];
221 for(Int_t i=0;i<3;i++){
222 for(Int_t j=i;j<3;j++){
223 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
227 fVertexXY->Fill(currentPos[0],currentPos[1]);
228 fVertexZ->Fill(currentPos[2]);
233 //______________________________________________________________________
234 Bool_t AliITSMeanVertexer::ComputeMean(){
235 // compute mean vertex
236 if(fNoEventsContr < 2) return kFALSE;
237 Double_t nevents = fNoEventsContr;
238 for(Int_t i=0;i<3;i++){
239 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
240 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
241 fAverPos[i] = fAverPosSum[i]/nevents;
243 for(Int_t i=0;i<3;i++){
244 for(Int_t j=i;j<3;j++){
245 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
246 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
250 fAverTracklets = fTotTracklets/nevents;
251 fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
252 fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
253 fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);