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 "AliITSRecPointContainer.h"
   9 #include "AliITSLoader.h"
  10 #include "AliLog.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"
  19 #include "AliMultiplicity.h"
  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(Bool_t mode):TObject(),
  43 fDetTypeRec(NULL),
  44 fVertexXY(NULL),
  45 fVertexZ(NULL),
  46 fNoEventsContr(0),
  47 fTotTracklets(0.),
  48 fAverTracklets(0.),
  49 fTotTrackletsSq(0.),
  50 fSigmaOnAverTracks(0.),
  51 fFilterOnContributors(0),
  52 fFilterOnTracklets(0),
  53 fMode(mode),
  54 fVertexer(NULL)
  55 {
  56   // Default Constructor
  57   for(Int_t i=0;i<3;i++){
  58     fWeighPosSum[i] = 0.;
  59     fWeighSigSum[i] = 0.;
  60     fAverPosSum[i] = 0.;
  61     fWeighPos[i] = 0.;
  62     fWeighSig[i] = 0.;
  63     fAverPos[i] = 0.;
  64     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  65     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
  66   }
  67
  68   // Histograms initialization
  69   const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
  70   const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
  71   fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
  72                        2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
  73                        2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
  74   fVertexXY->SetXTitle("X , cm");
  75   fVertexXY->SetYTitle("Y , cm");
  76   fVertexXY->SetOption("colz");
  77   fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
  78                        2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
  79   fVertexZ->SetXTitle("Z , cm");
  80 }
  81
  82 //______________________________________________________________________
  83 Bool_t AliITSMeanVertexer::Init() {
  84   // Initialize filters
  85   // Initialize geometry
  86   // Initialize ITS classes
  87
  88   AliGeomManager::LoadGeometry();
  89   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
  90
  91   AliITSInitGeometry initgeom;
  92   AliITSgeom *geom = initgeom.CreateAliITSgeom();
  93   if (!geom) return kFALSE;
  94   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
  95
  96   fDetTypeRec = new AliITSDetTypeRec();
  97   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
  98   fDetTypeRec->SetITSgeom(geom);
  99   fDetTypeRec->SetDefaults();
 100   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
 101
 102   // Initialize filter values to their defaults
 103   SetFilterOnContributors();
 104   SetFilterOnTracklets();
 105
 106   // Instatiate vertexer
 107   if (!fMode) {
 108     fVertexer = new AliITSVertexer3DTapan(1000);
 109   }
 110   else {
 111     fVertexer = new AliITSVertexer3D();
 112     fVertexer->SetDetTypeRec(fDetTypeRec);
 113     fVertexer->SetComputeMultiplicity(kTRUE);
 114   }
 115   return kTRUE;
 116 }
 117
 118 //______________________________________________________________________
 119 AliITSMeanVertexer::~AliITSMeanVertexer() {
 120   // Destructor
 121   delete fDetTypeRec;
 122   delete fVertexXY;
 123   delete fVertexZ;
 124   delete fVertexer;
 125 }
 126
 127 //______________________________________________________________________
 128 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
 129   // Performs SPD local reconstruction
 130   // and vertex finding
 131   // returns true in case a vertex is found
 132
 133   // Run SPD cluster finder
 134   AliITSRecPointContainer::Instance()->PrepareToRead();
 135   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
 136   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
 137
 138   Bool_t vtxOK = kFALSE;
 139   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
 140   if (!fMode) {
 141     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
 142   }
 143   else {
 144     AliMultiplicity *mult = fVertexer->GetMultiplicity();
 145     if(Filter(vtx,mult)) vtxOK = kTRUE;
 146   }
 147   delete clustersTree;
 148   if (vtxOK) AddToMean(vtx);
 149   if (vtx) delete vtx;
 150
 151   return vtxOK;
 152 }
 153
 154 //______________________________________________________________________
 155 void AliITSMeanVertexer::WriteVertices(const char *filename){
 156   // Compute mean vertex and
 157   // store it along with the histograms
 158   // in a file
 159
 160   TFile fmv(filename,"update");
 161
 162   if(ComputeMean()){
 163     Double_t cov[6];
 164     cov[0] =  fAverPosSq[0][0];  // variance x
 165     cov[1] =  fAverPosSq[0][1];  // cov xy
 166     cov[2] =  fAverPosSq[1][1];  // variance y
 167     cov[3] =  fAverPosSq[0][2];  // cov xz
 168     cov[4] =  fAverPosSq[1][2];  // cov yz
 169     cov[5] =  fAverPosSq[2][2];  // variance z
 170     AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
 171     mv.SetTitle("Mean Vertex");
 172     mv.SetName("MeanVertex");
 173     AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
 174     AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
 175     // we have to add chi2 here
 176     AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
 177
 178     mv.Write(mv.GetName(),TObject::kOverwrite);
 179     vtx.Write(vtx.GetName(),TObject::kOverwrite);
 180   }
 181   else {
 182     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
 183   }
 184
 185   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
 186   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
 187   fmv.Close();
 188 }
 189
 190 //______________________________________________________________________
 191 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
 192   // Apply selection criteria to events
 193   Bool_t status = kFALSE;
 194   if(!vert || !mult)return status;
 195   // Remove vertices reconstructed with vertexerZ
 196   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
 197   Int_t ncontr = vert->GetNContributors();
 198   Int_t ntracklets = mult->GetNumberOfTracklets();
 199   AliDebug(1,Form("Number of contributors = %d",ncontr));
 200   AliDebug(1,Form("Number of tracklets = %d",ntracklets));
 201   if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
 202   fTotTracklets += ntracklets;
 203   fTotTrackletsSq += ntracklets*ntracklets;
 204   return status;
 205 }
 206
 207 //______________________________________________________________________
 208 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
 209   // update mean vertex
 210   Double_t currentPos[3],currentSigma[3];
 211   vert->GetXYZ(currentPos);
 212   vert->GetSigmaXYZ(currentSigma);
 213   Bool_t goon = kTRUE;
 214   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
 215   if(!goon)return;
 216   for(Int_t i=0;i<3;i++){
 217     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
 218     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
 219     fAverPosSum[i]+=currentPos[i];
 220   }
 221   for(Int_t i=0;i<3;i++){
 222     for(Int_t j=i;j<3;j++){
 223       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
 224     }
 225   }
 226
 227   fVertexXY->Fill(currentPos[0],currentPos[1]);
 228   fVertexZ->Fill(currentPos[2]);
 229
 230   fNoEventsContr++;
 231 }
 232
 233 //______________________________________________________________________
 234 Bool_t AliITSMeanVertexer::ComputeMean(){
 235   // compute mean vertex
 236   if(fNoEventsContr < 2) return kFALSE;
 237   Double_t nevents = fNoEventsContr;
 238   for(Int_t i=0;i<3;i++){
 239     fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
 240     fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
 241     fAverPos[i] = fAverPosSum[i]/nevents;
 242   }
 243   for(Int_t i=0;i<3;i++){
 244     for(Int_t j=i;j<3;j++){
 245       fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
 246       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
 247     }
 248   }
 249
 250   fAverTracklets = fTotTracklets/nevents;
 251   fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
 252   fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
 253   fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
 254   return kTRUE;
 255 }