ITS/AliITSMeanVertexer.cxx

   1 #include <TFile.h>
   2 #include <TH1.h>
   3 #include <TH2.h>
   4 #include "AliGeomManager.h"
   5 #include "AliITSDetTypeRec.h"
   6 #include "AliITSInitGeometry.h"
   7 #include "AliITSMeanVertexer.h"
   8 #include "AliITSLoader.h"
   9 #include "AliLog.h"
  10 #include "AliRawReader.h"
  11 #include "AliRawReaderDate.h"
  12 #include "AliRawReaderRoot.h"
  13 #include "AliRunLoader.h"
  14 #include "AliITSVertexer3D.h"
  15 #include "AliITSVertexer3DTapan.h"
  16 #include "AliESDVertex.h"
  17 #include "AliMeanVertex.h"
  18 #include "AliMultiplicity.h"
  19
  20 ClassImp(AliITSMeanVertexer)
  21
  22 ///////////////////////////////////////////////////////////////////////
  23 //                                                                   //
  24 // Class to compute vertex position using SPD local reconstruction   //
  25 // An average vertex position using all the events                   //
  26 // is built and saved                                                //
  27 // Individual vertices can be optionally stored                      //
  28 // Origin: M.Masera  (masera@to.infn.it)                             //
  29 // Usage:
  30 // AliITSMeanVertexer mv("RawDataFileName");
  31 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
  32 // ....  optional setters ....
  33 // mv.Reconstruct();  // SPD local reconstruction
  34 // mv.DoVertices();
  35 // Resulting AliMeanVertex object is stored in a file named fMVFileName
  36 ///////////////////////////////////////////////////////////////////////
  37
  38 /* $Id$ */
  39
  40 //______________________________________________________________________
  41 AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
  42 fDetTypeRec(NULL),
  43 fVertexXY(NULL),
  44 fVertexZ(NULL),
  45 fNoEventsContr(0),
  46 fTotTracklets(0.),
  47 fAverTracklets(0.),
  48 fTotTrackletsSq(0.),
  49 fSigmaOnAverTracks(0.),
  50 fFilterOnContributors(0),
  51 fFilterOnTracklets(0)
  52 {
  53   // Default Constructor
  54   for(Int_t i=0;i<3;i++){
  55     fWeighPosSum[i] = 0.;
  56     fWeighSigSum[i] = 0.;
  57     fAverPosSum[i] = 0.;
  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     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
  63   }
  64
  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 }
  78
  79 //______________________________________________________________________
  80 Bool_t AliITSMeanVertexer::Init() {
  81   // Initialize filters
  82   // Initialize geometry
  83   // Initialize ITS classes
  84
  85   AliGeomManager::LoadGeometry();
  86   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
  87
  88   AliITSInitGeometry initgeom;
  89   AliITSgeom *geom = initgeom.CreateAliITSgeom();
  90   if (!geom) return kFALSE;
  91   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
  92
  93   fDetTypeRec = new AliITSDetTypeRec();
  94   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
  95   fDetTypeRec->SetITSgeom(geom);
  96   fDetTypeRec->SetDefaults();
  97   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
  98
  99   // Initialize filter values to their defaults
 100   SetFilterOnContributors();
 101   SetFilterOnTracklets();
 102
 103   return kTRUE;
 104 }
 105
 106 //______________________________________________________________________
 107 AliITSMeanVertexer::~AliITSMeanVertexer() {
 108   // Destructor
 109   delete fDetTypeRec;
 110   delete fVertexXY;
 111   delete fVertexZ;
 112 }
 113
 114 //______________________________________________________________________
 115 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
 116   // Performs SPD local reconstruction
 117   // and vertex finding
 118   // returns true in case a vertex is found
 119
 120   // Run SPD cluster finder
 121   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
 122   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
 123
 124   Bool_t vtxOK = kFALSE;
 125   AliESDVertex *vtx = NULL;
 126   // Run Tapan's vertex-finder
 127   if (!mode) {
 128     AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
 129     vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
 130     delete vertexer1;
 131     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
 132   }
 133   else {
 134   // Run standard vertexer3d
 135     AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
 136     vertexer2->SetDetTypeRec(fDetTypeRec);
 137     vertexer2->SetComputeMultiplicity(kTRUE);
 138     vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
 139     AliMultiplicity *mult = vertexer2->GetMultiplicity();
 140     delete vertexer2;
 141     if(Filter(vtx,mult)) vtxOK = kTRUE;
 142   }
 143   delete clustersTree;
 144   if (vtxOK) AddToMean(vtx);
 145   if (vtx) delete vtx;
 146
 147   return vtxOK;
 148 }
 149
 150 //______________________________________________________________________
 151 void AliITSMeanVertexer::WriteVertices(const char *filename){
 152   // Compute mean vertex and
 153   // store it along with the histograms
 154   // in a file
 155
 156   TFile fmv(filename,"update");
 157
 158   if(ComputeMean()){
 159     Double_t cov[6];
 160     cov[0] =  fAverPosSq[0][0];  // variance x
 161     cov[1] =  fAverPosSq[0][1];  // cov xy
 162     cov[2] =  fAverPosSq[1][1];  // variance y
 163     cov[3] =  fAverPosSq[0][2];  // cov xz
 164     cov[4] =  fAverPosSq[1][2];  // cov yz
 165     cov[5] =  fAverPosSq[2][2];  // variance z
 166     AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
 167     mv.SetTitle("Mean Vertex");
 168     mv.SetName("MeanVertex");
 169     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
 170     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
 171     // we have to add chi2 here
 172     AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
 173
 174     mv.Write(mv.GetName(),TObject::kOverwrite);
 175     vtx.Write(vtx.GetName(),TObject::kOverwrite);
 176   }
 177   else {
 178     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
 179   }
 180
 181   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
 182   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
 183   fmv.Close();
 184 }
 185
 186 //______________________________________________________________________
 187 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
 188   // Apply selection criteria to events
 189   Bool_t status = kFALSE;
 190   if(!vert || !mult)return status;
 191   // Remove vertices reconstructed with vertexerZ
 192   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
 193   Int_t ncontr = vert->GetNContributors();
 194   Int_t ntracklets = mult->GetNumberOfTracklets();
 195   AliDebug(1,Form("Number of contributors = %d",ncontr));
 196   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
 197   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
 198   fTotTracklets += ntracklets;
 199   fTotTrackletsSq += ntracklets*ntracklets;
 200   return status;
 201 }
 202
 203 //______________________________________________________________________
 204 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
 205   // update mean vertex
 206   Double_t currentPos[3],currentSigma[3];
 207   vert->GetXYZ(currentPos);
 208   vert->GetSigmaXYZ(currentSigma);
 209   Bool_t goon = kTRUE;
 210   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
 211   if(!goon)return;
 212   for(Int_t i=0;i<3;i++){
 213     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
 214     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
 215     fAverPosSum[i]+=currentPos[i];
 216   }
 217   for(Int_t i=0;i<3;i++){
 218     for(Int_t j=i;j<3;j++){
 219       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
 220     }
 221   }
 222
 223   fVertexXY->Fill(currentPos[0],currentPos[1]);
 224   fVertexZ->Fill(currentPos[2]);
 225
 226   fNoEventsContr++;
 227 }
 228
 229 //______________________________________________________________________
 230 Bool_t AliITSMeanVertexer::ComputeMean(){
 231   // compute mean vertex
 232   if(fNoEventsContr < 2) return kFALSE;
 233   Double_t nevents = fNoEventsContr;
 234   for(Int_t i=0;i<3;i++){
 235     fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
 236     fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
 237     fAverPos[i] = fAverPosSum[i]/nevents;
 238   }
 239   for(Int_t i=0;i<3;i++){
 240     for(Int_t j=i;j<3;j++){
 241       fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
 242       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
 243     }
 244   }
 245
 246   fAverTracklets = fTotTracklets/nevents;
 247   fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
 248   fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
 249   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
 250   return kTRUE;
 251 }