ITS/AliITSMeanVertexer.cxx

   1 #include <TFile.h>
   2 #include <TH1.h>
   3 #include <TH2.h>
   4 #include <TTree.h>
   5 #include "AliGeomManager.h"
   6 #include "AliITSDetTypeRec.h"
   7 #include "AliITSInitGeometry.h"
   8 #include "AliITSMeanVertexer.h"
   9 #include "AliITSRecPointContainer.h"
  10 #include "AliITSLoader.h"
  11 #include "AliLog.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"
  20
  21 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
  22 ClassImp(AliITSMeanVertexer)
  23
  24 ///////////////////////////////////////////////////////////////////////
  25 //                                                                   //
  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)                             //
  31 // Usage:
  32 // AliITSMeanVertexer mv("RawDataFileName");
  33 // mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
  34 // ....  optional setters ....
  35 // mv.Reconstruct();  // SPD local reconstruction
  36 // mv.DoVertices();
  37 // Resulting AliMeanVertex object is stored in a file named fMVFileName
  38 ///////////////////////////////////////////////////////////////////////
  39
  40 /* $Id$ */
  41
  42 //______________________________________________________________________
  43 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
  44 fDetTypeRec(NULL),
  45 fVertexXY(NULL),
  46 fVertexZ(NULL),
  47 fNoEventsContr(0),
  48 fTotContributors(0.),
  49 fAverContributors(0.),
  50 fFilterOnContributors(0),
  51 fMode(mode),
  52 fVertexer(NULL),
  53 fAccEvents(fgkMaxNumOfEvents),
  54 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
  55 fClu0(NULL),
  56 fIndex(0),
  57 fErrXCut(0.),
  58 fRCut(0.),
  59 fLowSPD0(0),
  60 fHighSPD0(0)
  61 {
  62   // Default Constructor
  63   Reset();
  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   fClu0 = new UInt_t [fgkMaxNumOfEvents];
  78   for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
  79   fAccEvents.ResetAllBits(kTRUE);
  80 }
  81
  82 //______________________________________________________________________
  83 void AliITSMeanVertexer::Reset(){
  84   // Reset averages
  85   for(Int_t i=0;i<3;i++){
  86     fWeighPosSum[i] = 0.;
  87     fWeighSigSum[i] = 0.;
  88     fAverPosSum[i] = 0.;
  89     fWeighPos[i] = 0.;
  90     fWeighSig[i] = 0.;
  91     fAverPos[i] = 0.;
  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.;
  94     fNoEventsContr=0;
  95     fTotContributors = 0.;
  96     if(fVertexXY)fVertexXY->Reset();
  97     if(fVertexZ)fVertexZ->Reset();
  98   }
  99 }
 100 //______________________________________________________________________
 101 Bool_t AliITSMeanVertexer::Init() {
 102   // Initialize filters
 103   // Initialize geometry
 104   // Initialize ITS classes
 105
 106   AliGeomManager::LoadGeometry();
 107   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
 108
 109   AliITSInitGeometry initgeom;
 110   AliITSgeom *geom = initgeom.CreateAliITSgeom();
 111   if (!geom) return kFALSE;
 112   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
 113
 114   fDetTypeRec = new AliITSDetTypeRec();
 115   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
 116   fDetTypeRec->SetITSgeom(geom);
 117   fDetTypeRec->SetDefaults();
 118   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
 119
 120   // Initialize filter values to their defaults
 121   SetFilterOnContributors();
 122   SetCutOnErrX();
 123   SetCutOnR();
 124   SetCutOnCls();
 125
 126   // Instatiate vertexer
 127   if (!fMode) {
 128     fVertexer = new AliITSVertexer3DTapan(1000);
 129   }
 130   else {
 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);
 140   }
 141   return kTRUE;
 142 }
 143
 144 //______________________________________________________________________
 145 AliITSMeanVertexer::~AliITSMeanVertexer() {
 146   // Destructor
 147   delete fDetTypeRec;
 148   delete fVertexXY;
 149   delete fVertexZ;
 150   delete fVertexer;
 151 }
 152
 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
 158
 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");
 163
 164   Bool_t vtxOK = kFALSE;
 165   UInt_t nl1=0;
 166   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
 167   if (!fMode) {
 168     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
 169   }
 170   else {
 171     AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
 172     nl1=rpcont->GetNClustersInLayerFast(1);
 173     if(Filter(vtx,nl1)) vtxOK = kTRUE;
 174   }
 175   delete clustersTree;
 176   if (vtxOK && fMode){
 177     new(fVertArray[fIndex])AliESDVertex(*vtx);
 178     fClu0[fIndex]=nl1;
 179     fIndex++;
 180   }
 181   if (vtx) delete vtx;
 182
 183   return vtxOK;
 184 }
 185
 186 //______________________________________________________________________
 187 void AliITSMeanVertexer::WriteVertices(const char *filename){
 188   // Compute mean vertex and
 189   // store it along with the histograms
 190   // in a file
 191
 192   TFile fmv(filename,"update");
 193   Bool_t itisOK = kFALSE;
 194   if(ComputeMean(kFALSE)){
 195     if(ComputeMean(kTRUE)){
 196       Double_t cov[6];
 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");
 208
 209       mv.Write(mv.GetName(),TObject::kOverwrite);
 210       vtx.Write(vtx.GetName(),TObject::kOverwrite);
 211       itisOK = kTRUE;
 212     }
 213   }
 214
 215   if(!itisOK){
 216     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
 217   }
 218
 219   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
 220   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
 221   fmv.Close();
 222 }
 223
 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)) {
 236     status = kTRUE;
 237     fAccEvents.SetBitNumber(fIndex);
 238   }
 239   return status;
 240 }
 241
 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);
 248   Bool_t goon = kTRUE;
 249   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
 250   if(!goon)return;
 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];
 255   }
 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];
 259     }
 260   }
 261
 262   fTotContributors+=vert->GetNContributors();
 263   fVertexXY->Fill(currentPos[0],currentPos[1]);
 264   fVertexZ->Fill(currentPos[2]);
 265
 266   fNoEventsContr++;
 267 }
 268
 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
 274    Double_t wpos[3];
 275    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
 276    Reset();
 277    for(Int_t i=0;i<fVertArray.GetEntries();i++){
 278      AliESDVertex* vert = NULL;
 279
 280      // second pass
 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);
 287      }
 288
 289      if(!fAccEvents.TestBitNumber(i))continue;
 290      vert=(AliESDVertex*)fVertArray[i];
 291      AddToMean(vert);
 292    }
 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;
 299    }
 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];
 304      }
 305    }
 306
 307    fAverContributors = fTotContributors/nevents;
 308    return kTRUE;
 309  }