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"
20 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
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(),
48 fAverContributors(0.),
49 fFilterOnContributors(0),
52 fAccEvents(fgkMaxNumOfEvents),
53 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
61 // Default Constructor
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);
81 //______________________________________________________________________
82 void AliITSMeanVertexer::Reset(){
84 for(Int_t i=0;i<3;i++){
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.;
94 fTotContributors = 0.;
95 if(fVertexXY)fVertexXY->Reset();
96 if(fVertexZ)fVertexZ->Reset();
99 //______________________________________________________________________
100 Bool_t AliITSMeanVertexer::Init() {
101 // Initialize filters
102 // Initialize geometry
103 // Initialize ITS classes
105 AliGeomManager::LoadGeometry();
106 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
108 AliITSInitGeometry initgeom;
109 AliITSgeom *geom = initgeom.CreateAliITSgeom();
110 if (!geom) return kFALSE;
111 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
113 fDetTypeRec = new AliITSDetTypeRec();
114 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
115 fDetTypeRec->SetITSgeom(geom);
116 fDetTypeRec->SetDefaults();
117 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
119 // Initialize filter values to their defaults
120 SetFilterOnContributors();
125 // Instatiate vertexer
127 fVertexer = new AliITSVertexer3DTapan(1000);
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);
143 //______________________________________________________________________
144 AliITSMeanVertexer::~AliITSMeanVertexer() {
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
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");
163 Bool_t vtxOK = kFALSE;
165 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
167 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
170 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
171 nl1=rpcont->GetNClustersInLayerFast(1);
172 if(Filter(vtx,nl1)) vtxOK = kTRUE;
176 new(fVertArray[fIndex])AliESDVertex(*vtx);
185 //______________________________________________________________________
186 void AliITSMeanVertexer::WriteVertices(const char *filename){
187 // Compute mean vertex and
188 // store it along with the histograms
191 TFile fmv(filename,"update");
192 Bool_t itisOK = kFALSE;
193 if(ComputeMean(kFALSE)){
194 if(ComputeMean(kTRUE)){
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");
208 mv.Write(mv.GetName(),TObject::kOverwrite);
209 vtx.Write(vtx.GetName(),TObject::kOverwrite);
215 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
218 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
219 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
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)) {
236 fAccEvents.SetBitNumber(fIndex);
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);
248 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
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];
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];
261 fTotContributors+=vert->GetNContributors();
262 fVertexXY->Fill(currentPos[0],currentPos[1]);
263 fVertexZ->Fill(currentPos[2]);
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
274 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
276 for(Int_t i=0;i<fVertArray.GetEntries();i++){
277 AliESDVertex* vert = NULL;
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);
288 if(!fAccEvents.TestBitNumber(i))continue;
289 vert=(AliESDVertex*)fVertArray[i];
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;
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];
306 fAverContributors = fTotContributors/nevents;