5 #include "AliGeomManager.h"
6 #include "AliITSDetTypeRec.h"
7 #include "AliITSInitGeometry.h"
8 #include "AliITSMeanVertexer.h"
9 #include "AliITSRecPointContainer.h"
10 #include "AliITSLoader.h"
12 #include "AliRawReader.h"
13 #include "AliRawReaderDate.h"
14 #include "AliRawReaderRoot.h"
15 #include "AliRunLoader.h"
16 #include "AliITSVertexer3D.h"
17 #include "AliITSVertexer3DTapan.h"
18 #include "AliESDVertex.h"
19 #include "AliMeanVertex.h"
21 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
22 ClassImp(AliITSMeanVertexer)
24 ///////////////////////////////////////////////////////////////////////
26 // Class to compute vertex position using SPD local reconstruction //
27 // An average vertex position using all the events //
28 // is built and saved //
29 // Individual vertices can be optionally stored //
30 // Origin: M.Masera (masera@to.infn.it) //
32 // AliITSMeanVertexer mv("RawDataFileName");
33 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
34 // .... optional setters ....
35 // mv.Reconstruct(); // SPD local reconstruction
37 // Resulting AliMeanVertex object is stored in a file named fMVFileName
38 ///////////////////////////////////////////////////////////////////////
42 //______________________________________________________________________
43 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
49 fAverContributors(0.),
50 fFilterOnContributors(0),
53 fAccEvents(fgkMaxNumOfEvents),
54 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
62 // Default Constructor
65 // Histograms initialization
66 const Float_t xLimit = 5.0, yLimit = 5.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 fVertexXY->SetXTitle("X , cm");
72 fVertexXY->SetYTitle("Y , cm");
73 fVertexXY->SetOption("colz");
74 fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
75 2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
76 fVertexZ->SetXTitle("Z , cm");
77 fClu0 = new UInt_t [fgkMaxNumOfEvents];
78 for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
79 fAccEvents.ResetAllBits(kTRUE);
82 //______________________________________________________________________
83 void AliITSMeanVertexer::Reset(){
85 for(Int_t i=0;i<3;i++){
92 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
93 for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
95 fTotContributors = 0.;
96 if(fVertexXY)fVertexXY->Reset();
97 if(fVertexZ)fVertexZ->Reset();
100 //______________________________________________________________________
101 Bool_t AliITSMeanVertexer::Init() {
102 // Initialize filters
103 // Initialize geometry
104 // Initialize ITS classes
106 AliGeomManager::LoadGeometry();
107 if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
109 AliITSInitGeometry initgeom;
110 AliITSgeom *geom = initgeom.CreateAliITSgeom();
111 if (!geom) return kFALSE;
112 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
114 fDetTypeRec = new AliITSDetTypeRec();
115 fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
116 fDetTypeRec->SetITSgeom(geom);
117 fDetTypeRec->SetDefaults();
118 fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
120 // Initialize filter values to their defaults
121 SetFilterOnContributors();
126 // Instatiate vertexer
128 fVertexer = new AliITSVertexer3DTapan(1000);
131 fVertexer = new AliITSVertexer3D();
132 fVertexer->SetDetTypeRec(fDetTypeRec);
133 AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
134 alias->SetWideFiducialRegion(40.,0.5);
135 alias->SetNarrowFiducialRegion(0.5,0.5);
136 alias->SetDeltaPhiCuts(0.5,0.025);
137 alias->SetDCACut(0.1);
138 alias->SetPileupAlgo(3);
139 fVertexer->SetComputeMultiplicity(kFALSE);
144 //______________________________________________________________________
145 AliITSMeanVertexer::~AliITSMeanVertexer() {
153 //______________________________________________________________________
154 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
155 // Performs SPD local reconstruction
156 // and vertex finding
157 // returns true in case a vertex is found
159 // Run SPD cluster finder
160 AliITSRecPointContainer::Instance()->PrepareToRead();
161 TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
162 fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
164 Bool_t vtxOK = kFALSE;
166 AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
168 if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
171 AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
172 nl1=rpcont->GetNClustersInLayerFast(1);
173 if(Filter(vtx,nl1)) vtxOK = kTRUE;
177 new(fVertArray[fIndex])AliESDVertex(*vtx);
186 //______________________________________________________________________
187 void AliITSMeanVertexer::WriteVertices(const char *filename){
188 // Compute mean vertex and
189 // store it along with the histograms
192 TFile fmv(filename,"update");
193 Bool_t itisOK = kFALSE;
194 if(ComputeMean(kFALSE)){
195 if(ComputeMean(kTRUE)){
197 cov[0] = fAverPosSq[0][0]; // variance x
198 cov[1] = fAverPosSq[0][1]; // cov xy
199 cov[2] = fAverPosSq[1][1]; // variance y
200 cov[3] = fAverPosSq[0][2]; // cov xz
201 cov[4] = fAverPosSq[1][2]; // cov yz
202 cov[5] = fAverPosSq[2][2]; // variance z
203 AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
204 mv.SetTitle("Mean Vertex");
205 mv.SetName("MeanVertex");
206 // we have to add chi2 here
207 AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
209 mv.Write(mv.GetName(),TObject::kOverwrite);
210 vtx.Write(vtx.GetName(),TObject::kOverwrite);
216 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
219 fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
220 fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
224 //______________________________________________________________________
225 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
226 // Apply selection criteria to events
227 Bool_t status = kFALSE;
228 if(!vert)return status;
229 // Remove vertices reconstructed with vertexerZ
230 if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
231 Int_t ncontr = vert->GetNContributors();
232 AliDebug(1,Form("Number of contributors = %d",ncontr));
233 Double_t ex = vert->GetXRes();
234 if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
235 ((ex*1000.)<fErrXCut)) {
237 fAccEvents.SetBitNumber(fIndex);
242 //______________________________________________________________________
243 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
244 // update mean vertex
245 Double_t currentPos[3],currentSigma[3];
246 vert->GetXYZ(currentPos);
247 vert->GetSigmaXYZ(currentSigma);
249 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
251 for(Int_t i=0;i<3;i++){
252 fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
253 fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
254 fAverPosSum[i]+=currentPos[i];
256 for(Int_t i=0;i<3;i++){
257 for(Int_t j=i;j<3;j++){
258 fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
262 fTotContributors+=vert->GetNContributors();
263 fVertexXY->Fill(currentPos[0],currentPos[1]);
264 fVertexZ->Fill(currentPos[2]);
269 //______________________________________________________________________
270 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
271 // compute mean vertex
272 // called once with killOutliers = kFALSE and then again with
273 // killOutliers = kTRUE
275 for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
277 for(Int_t i=0;i<fVertArray.GetEntries();i++){
278 AliESDVertex* vert = NULL;
281 if(killOutliers && fAccEvents.TestBitNumber(i)){
282 vert=(AliESDVertex*)fVertArray[i];
283 Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
284 dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
285 dist=sqrt(dist)*10.; // distance in mm
286 if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
289 if(!fAccEvents.TestBitNumber(i))continue;
290 vert=(AliESDVertex*)fVertArray[i];
293 if(fNoEventsContr < 2) return kFALSE;
294 Double_t nevents = fNoEventsContr;
295 for(Int_t i=0;i<3;i++){
296 fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
297 fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
298 fAverPos[i] = fAverPosSum[i]/nevents;
300 for(Int_t i=0;i<3;i++){
301 for(Int_t j=i;j<3;j++){
302 fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
303 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
307 fAverContributors = fTotContributors/nevents;