ITS/AliITSMeanVertexer.cxx

   1 #include <TFile.h>
   2 #include "AliGeomManager.h"
   3 #include "AliHeader.h"
   4 #include "AliITSDetTypeRec.h"
   5 #include "AliITSInitGeometry.h"
   6 #include "AliITSMeanVertexer.h"
   7 #include "AliITSLoader.h"
   8 #include "AliLog.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"
  17
  18 const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root";
  19
  20
  21 ClassImp(AliITSMeanVertexer)
  22
  23 ///////////////////////////////////////////////////////////////////////
  24 //                                                                   //
  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)                             //
  30 // Usage:
  31 // AliITSMeanVertexer mv("RawDataFileName");
  32 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
  33 // ....  optional setters ....
  34 // mv.Reconstruct();  // SPD local reconstruction
  35 // mv.DoVertices();
  36 // Resulting AliMeanVertex object is stored in a file named fMVFileName
  37 ///////////////////////////////////////////////////////////////////////
  38
  39 /* $Id$ */
  40
  41 //______________________________________________________________________
  42 AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
  43 fLoaderFileName(),
  44 fGeometryFileName(),
  45 fMVFileName(),
  46 fRawReader(),
  47 fRunLoader(),
  48 fNoEventsContr(0),
  49 fTotTracklets(0.),
  50 fAverTracklets(0.),
  51 fSigmaOnAverTracks(0.),
  52 fFilterOnContributors(0),
  53 fFilterOnTracklets(0),
  54 fWriteVertices(kFALSE)
  55 {
  56   // Default Constructor
  57   for(Int_t i=0;i<3;i++){
  58     fWeighPos[i] = 0.;
  59     fWeighSig[i] = 0.;
  60     fAverPos[i] = 0.;
  61     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  62   }
  63
  64 }
  65
  66 //______________________________________________________________________
  67 AliITSMeanVertexer::AliITSMeanVertexer(TString filename):TObject(),
  68 fLoaderFileName(),
  69 fGeometryFileName(),
  70 fMVFileName(fgkMVFileNameDefault),
  71 fRawReader(),
  72 fRunLoader(),
  73 fNoEventsContr(0),
  74 fTotTracklets(0.),
  75 fAverTracklets(0.),
  76 fSigmaOnAverTracks(0.),
  77 fFilterOnContributors(0),
  78 fFilterOnTracklets(0),
  79 fWriteVertices(kTRUE)
  80 {
  81   // Standard constructor
  82
  83   for(Int_t i=0;i<3;i++){
  84     fWeighPos[i] = 0.;
  85     fWeighSig[i] = 0.;
  86     fAverPos[i] = 0.;
  87     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  88   }
  89   SetLoaderFileName();
  90   SetGeometryFileName();
  91
  92   Init(filename);
  93 }
  94 //______________________________________________________________________
  95 AliITSMeanVertexer::AliITSMeanVertexer(TString filename,
  96                                        TString loaderfilename,
  97                                        TString geometryfilename):TObject(),
  98 fLoaderFileName(),
  99 fGeometryFileName(),
 100 fMVFileName(fgkMVFileNameDefault),
 101 fRawReader(),
 102 fRunLoader(),
 103 fNoEventsContr(0),
 104 fTotTracklets(0.),
 105 fAverTracklets(0.),
 106 fSigmaOnAverTracks(0.),
 107 fFilterOnContributors(0),
 108 fFilterOnTracklets(0),
 109 fWriteVertices(kTRUE)
 110 {
 111   // Standard constructor with explicit geometry file name assignment
 112   for(Int_t i=0;i<3;i++){
 113     fWeighPos[i] = 0.;
 114     fWeighSig[i] = 0.;
 115     fAverPos[i] = 0.;
 116     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
 117   }
 118   SetLoaderFileName(loaderfilename);
 119   SetGeometryFileName(geometryfilename);
 120   Init(filename);
 121 }
 122
 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");
 128   }
 129   // if file name ends with root a raw reader ROOT is assumed
 130   if(filename.EndsWith(".root")){
 131     fRawReader = new AliRawReaderRoot(filename);
 132   }
 133   else {  // DATE raw reader is assumed
 134     fRawReader = new AliRawReaderDate(filename);
 135     fRawReader->SelectEvents(7);
 136   }
 137   fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
 138   fRunLoader->MakeTree("E");
 139   Int_t iEvent = 0;
 140   while (fRawReader->NextEvent()) {
 141     fRunLoader->SetEventNumber(iEvent);
 142     fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 143                                    iEvent, iEvent);
 144     fRunLoader->MakeTree("H");
 145     fRunLoader->TreeE()->Fill();
 146     iEvent++;
 147   }
 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
 158
 159   AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 160
 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();
 168 }
 169
 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)
 184 {
 185   // Copy constructor
 186   // Copies are not allowed. The method is protected to avoid misuse.
 187   AliFatal("Copy constructor not allowed\n");
 188
 189 }
 190
 191 //______________________________________________________________________
 192 AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer&  /* vtxr */ ){
 193   // Assignment operator
 194   AliError("Assignment operator not allowed\n");
 195   return *this;
 196 }
 197
 198 //______________________________________________________________________
 199 AliITSMeanVertexer::~AliITSMeanVertexer() {
 200   // Destructor
 201   delete fRawReader;
 202   delete fRunLoader;
 203
 204 }
 205
 206 //______________________________________________________________________
 207 void AliITSMeanVertexer::Reconstruct(){
 208   // Performs SPD local reconstruction
 209   AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
 210   if (!loader) {
 211     AliFatal("ITS loader not found");
 212     return;
 213   }
 214   AliITSDetTypeRec* rec = new AliITSDetTypeRec();
 215   rec->SetITSgeom(loader->GetITSgeom());
 216   rec->SetDefaults();
 217
 218   rec->SetDefaultClusterFindersV2(kTRUE);
 219
 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();
 230     if(!tR){
 231       AliFatal("Tree R pointer not found - Abort \n");
 232       break;
 233     }
 234     rec->DigitsToRecPoints(fRawReader,tR,"SPD");
 235     rec->ResetRecPoints();
 236     rec->ResetClusters();
 237     loader->WriteRecPoints("OVERWRITE");
 238     loader->UnloadRecPoints();
 239   }
 240
 241 }
 242
 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("ITS.Vert3D.root");
 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     vert = dovert->FindVertexForCurrentEvent(i);
 253     AliMultiplicity *mult = dovert->GetMultiplicity();
 254     if(Filter(vert,mult)){
 255       AddToMean(vert);
 256       if(fWriteVertices){
 257         loader->PostVertex(vert);
 258         loader->WriteVertices();
 259       }
 260     }
 261     else {
 262       if(vert)delete vert;
 263     }
 264   }
 265   if(ComputeMean()){
 266     Double_t cov[6];
 267     cov[0] =  fAverPosSq[0][0];  // variance x
 268     cov[1] =  fAverPosSq[0][1];  // cov xy
 269     cov[2] =  fAverPosSq[1][1];  // variance y
 270     cov[3] =  fAverPosSq[0][2];  // cov xz
 271     cov[4] =  fAverPosSq[1][2];  // cov yz
 272     cov[5] =  fAverPosSq[2][2];  // variance z
 273     AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
 274     mv.SetTitle("Mean Vertex");
 275     mv.SetName("Meanvertex");
 276     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
 277     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
 278     TFile fmv(fMVFileName.Data(),"recreate");
 279     mv.Write();
 280     fmv.Close();
 281   }
 282   else {
 283     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
 284   }
 285   delete dovert;
 286 }
 287
 288 //______________________________________________________________________
 289 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
 290   // Apply selection criteria to events
 291   Bool_t status = kFALSE;
 292   if(!vert || !mult)return status;
 293   Int_t ncontr = vert->GetNContributors();
 294   Int_t ntracklets = mult->GetNumberOfTracklets();
 295   AliDebug(1,Form("Number of contributors = %d",ncontr));
 296   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
 297   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
 298   fAverTracklets += ntracklets;
 299   fSigmaOnAverTracks += ntracklets*ntracklets;
 300   return status;
 301 }
 302
 303 //______________________________________________________________________
 304 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
 305   // update mean vertex
 306   Double_t currentPos[3],currentSigma[3];
 307   vert->GetXYZ(currentPos);
 308   vert->GetSigmaXYZ(currentSigma);
 309   Bool_t goon = kTRUE;
 310   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
 311   if(!goon)return;
 312   for(Int_t i=0;i<3;i++){
 313     fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
 314     fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
 315     fAverPos[i]+=currentPos[i];
 316   }
 317   for(Int_t i=0;i<3;i++){
 318     for(Int_t j=i;j<3;j++){
 319       fAverPosSq[i][j] += currentPos[i] * currentPos[j];
 320     }
 321   }
 322   fNoEventsContr++;
 323 }
 324
 325 //______________________________________________________________________
 326 Bool_t AliITSMeanVertexer::ComputeMean(){
 327   // compute mean vertex
 328   if(fNoEventsContr < 2) return kFALSE;
 329   Double_t nevents = fNoEventsContr;
 330   for(Int_t i=0;i<3;i++){
 331     fWeighPos[i] /= fWeighSig[i];
 332     fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
 333     fAverPos[i] /= nevents;
 334   }
 335   for(Int_t i=0;i<3;i++){
 336     for(Int_t j=i;j<3;j++){
 337       fAverPosSq[i][j] /= (nevents -1.);
 338       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
 339     }
 340   }
 341   fTotTracklets = fAverTracklets;  //  total number of tracklets used
 342   fAverTracklets /= nevents;
 343   fSigmaOnAverTracks /= (nevents - 1);
 344   Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
 345   fSigmaOnAverTracks -= tmp;
 346   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
 347   return kTRUE;
 348 }