2 #include "AliGeomManager.h"
4 #include "AliITSDetTypeRec.h"
5 #include "AliITSInitGeometry.h"
6 #include "AliITSMeanVertexer.h"
7 #include "AliITSLoader.h"
9 #include "AliRawReader.h"
10 #include "AliRawReaderDate.h"
11 #include "AliRawReaderRoot.h"
12 #include "AliRunLoader.h"
13 #include "AliITSVertexer3D.h"
14 #include "AliESDVertex.h"
15 #include "AliMeanVertex.h"
16 #include "AliMultiplicity.h"
18 const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root";
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():TObject(),
51 fSigmaOnAverTracks(0.),
52 fFilterOnContributors(0),
53 fFilterOnTracklets(0),
54 fWriteVertices(kFALSE)
56 // Default Constructor
57 for(Int_t i=0;i<3;i++){
61 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
66 //______________________________________________________________________
67 AliITSMeanVertexer::AliITSMeanVertexer(TString &filename):TObject(),
70 fMVFileName(fgkMVFileNameDefault),
76 fSigmaOnAverTracks(0.),
77 fFilterOnContributors(0),
78 fFilterOnTracklets(0),
81 // Standard constructor
83 for(Int_t i=0;i<3;i++){
87 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
90 SetGeometryFileName();
94 //______________________________________________________________________
95 AliITSMeanVertexer::AliITSMeanVertexer(TString &filename,
96 TString &loaderfilename,
97 TString &geometryfilename):TObject(),
100 fMVFileName(fgkMVFileNameDefault),
106 fSigmaOnAverTracks(0.),
107 fFilterOnContributors(0),
108 fFilterOnTracklets(0),
109 fWriteVertices(kTRUE)
111 // Standard constructor with explicit geometry file name assignment
112 for(Int_t i=0;i<3;i++){
116 for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
118 SetLoaderFileName(loaderfilename);
119 SetGeometryFileName(geometryfilename);
123 //______________________________________________________________________
124 void AliITSMeanVertexer::Init(TString &filename){
125 // Initialization part common to different constructors
126 if(filename.IsNull()){
127 AliFatal("Please, provide a valid file name for raw data file\n");
129 // if file name ends with root a raw reader ROOT is assumed
130 if(filename.EndsWith(".root")){
131 fRawReader = new AliRawReaderRoot(filename);
133 else { // DATE raw reader is assumed
134 fRawReader = new AliRawReaderDate(filename);
135 fRawReader->SelectEvents(7);
137 fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
138 fRunLoader->MakeTree("E");
140 while (fRawReader->NextEvent()) {
141 fRunLoader->SetEventNumber(iEvent);
142 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
144 fRunLoader->MakeTree("H");
145 fRunLoader->TreeE()->Fill();
148 fRawReader->RewindEvents();
149 fRunLoader->SetNumberOfEventsPerFile(iEvent);
150 fRunLoader->WriteHeader("OVERWRITE");
151 Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS");
152 if(retval != 0)AliFatal("Not able to add ITS detector");
153 AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName());
154 fRunLoader->AddLoader(loader);
155 fRunLoader->CdGAFile();
156 fRunLoader->Write(0, TObject::kOverwrite);
157 // Initialize geometry
159 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
161 AliITSInitGeometry initgeom;
162 AliITSgeom *geom = initgeom.CreateAliITSgeom();
163 printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
164 loader->SetITSgeom(geom);
165 // Initialize filter values to their defaults
166 SetFilterOnContributors();
167 SetFilterOnTracklets();
170 //______________________________________________________________________
171 AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr),
172 fLoaderFileName(vtxr.fLoaderFileName),
173 fGeometryFileName(vtxr.fGeometryFileName),
174 fMVFileName(vtxr.fMVFileName),
175 fRawReader(vtxr.fRawReader),
176 fRunLoader(vtxr.fRunLoader),
177 fNoEventsContr(vtxr.fNoEventsContr),
178 fTotTracklets(vtxr.fTotTracklets),
179 fAverTracklets(vtxr.fAverTracklets),
180 fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks),
181 fFilterOnContributors(vtxr.fFilterOnContributors),
182 fFilterOnTracklets(vtxr.fFilterOnTracklets),
183 fWriteVertices(vtxr.fWriteVertices)
186 // Copies are not allowed. The method is protected to avoid misuse.
187 AliFatal("Copy constructor not allowed\n");
191 //______________________________________________________________________
192 AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){
193 // Assignment operator
194 AliError("Assignment operator not allowed\n");
198 //______________________________________________________________________
199 AliITSMeanVertexer::~AliITSMeanVertexer() {
206 //______________________________________________________________________
207 void AliITSMeanVertexer::Reconstruct(){
208 // Performs SPD local reconstruction
209 AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
211 AliFatal("ITS loader not found");
214 AliITSDetTypeRec* rec = new AliITSDetTypeRec();
215 rec->SetITSgeom(loader->GetITSgeom());
218 rec->SetDefaultClusterFindersV2(kTRUE);
220 Int_t nEvents = fRunLoader->GetNumberOfEvents();
221 fRawReader->RewindEvents();
222 for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
223 fRawReader->NextEvent();
224 fRunLoader->GetEvent(iEvent);
225 AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent));
226 loader->LoadRecPoints("update");
227 loader->CleanRecPoints();
228 loader->MakeRecPointsContainer();
229 TTree *tR = loader->TreeR();
231 AliFatal("Tree R pointer not found - Abort \n");
234 rec->DigitsToRecPoints(fRawReader,tR,"SPD");
235 rec->ResetRecPoints();
236 rec->ResetClusters();
237 loader->WriteRecPoints("OVERWRITE");
238 loader->UnloadRecPoints();
243 //______________________________________________________________________
244 void AliITSMeanVertexer::DoVertices(){
245 // Loop on all events and compute 3D vertices
246 AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
247 AliITSVertexer3D *dovert = new AliITSVertexer3D();
248 AliESDVertex *vert = 0;
249 Int_t nevents = fRunLoader->TreeE()->GetEntries();
250 for(Int_t i=0; i<nevents; i++){
251 fRunLoader->GetEvent(i);
252 TTree* cltree = loader->TreeR();
253 vert = dovert->FindVertexForCurrentEvent(cltree);
254 AliMultiplicity *mult = dovert->GetMultiplicity();
255 if(Filter(vert,mult)){
258 loader->PostVertex(vert);
259 loader->WriteVertices();
268 cov[0] = fAverPosSq[0][0]; // variance x
269 cov[1] = fAverPosSq[0][1]; // cov xy
270 cov[2] = fAverPosSq[1][1]; // variance y
271 cov[3] = fAverPosSq[0][2]; // cov xz
272 cov[4] = fAverPosSq[1][2]; // cov yz
273 cov[5] = fAverPosSq[2][2]; // variance z
274 AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
275 mv.SetTitle("Mean Vertex");
276 mv.SetName("Meanvertex");
277 AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
278 AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
279 TFile fmv(fMVFileName.Data(),"recreate");
284 AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
289 //______________________________________________________________________
290 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
291 // Apply selection criteria to events
292 Bool_t status = kFALSE;
293 if(!vert || !mult)return status;
294 Int_t ncontr = vert->GetNContributors();
295 Int_t ntracklets = mult->GetNumberOfTracklets();
296 AliDebug(1,Form("Number of contributors = %d",ncontr));
297 AliDebug(1,Form("Number of tracklets = %d",ntracklets));
298 if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
299 fAverTracklets += ntracklets;
300 fSigmaOnAverTracks += ntracklets*ntracklets;
304 //______________________________________________________________________
305 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
306 // update mean vertex
307 Double_t currentPos[3],currentSigma[3];
308 vert->GetXYZ(currentPos);
309 vert->GetSigmaXYZ(currentSigma);
311 for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
313 for(Int_t i=0;i<3;i++){
314 fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
315 fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
316 fAverPos[i]+=currentPos[i];
318 for(Int_t i=0;i<3;i++){
319 for(Int_t j=i;j<3;j++){
320 fAverPosSq[i][j] += currentPos[i] * currentPos[j];
326 //______________________________________________________________________
327 Bool_t AliITSMeanVertexer::ComputeMean(){
328 // compute mean vertex
329 if(fNoEventsContr < 2) return kFALSE;
330 Double_t nevents = fNoEventsContr;
331 for(Int_t i=0;i<3;i++){
332 fWeighPos[i] /= fWeighSig[i];
333 fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
334 fAverPos[i] /= nevents;
336 for(Int_t i=0;i<3;i++){
337 for(Int_t j=i;j<3;j++){
338 fAverPosSq[i][j] /= (nevents -1.);
339 fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
342 fTotTracklets = fAverTracklets; // total number of tracklets used
343 fAverTracklets /= nevents;
344 fSigmaOnAverTracks /= (nevents - 1);
345 Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
346 fSigmaOnAverTracks -= tmp;
347 fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);