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 //#include "AliCodeTimer.h"
  21
  22 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
  23 ClassImp(AliITSMeanVertexer)
  24
  25 ///////////////////////////////////////////////////////////////////////
  26 //                                                                   //
  27 // Class to compute vertex position using SPD local reconstruction   //
  28 // An average vertex position using all the events                   //
  29 // is built and saved                                                //
  30 // Individual vertices can be optionally stored                      //
  31 // Origin: M.Masera  (masera@to.infn.it)                             //
  32 // Usage:                                                            //
  33 // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm    //
  34 ///////////////////////////////////////////////////////////////////////
  35
  36 /* $Id$ */
  37
  38 //______________________________________________________________________
  39 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
  40 fDetTypeRec(NULL),
  41 fVertexXY(NULL),
  42 fVertexZ(NULL),
  43 fNoEventsContr(0),
  44 fTotContributors(0.),
  45 fAverContributors(0.),
  46 fFilterOnContributors(0),
  47 fMode(mode),
  48 fVertexer(NULL),
  49 fAccEvents(fgkMaxNumOfEvents),
  50 fVertArray("AliESDVertex",fgkMaxNumOfEvents),
  51 fClu0(NULL),
  52 fIndex(0),
  53 fErrXCut(0.),
  54 fRCut(0.),
  55 fLowSPD0(0),
  56 fHighSPD0(0),
  57 fMultH(NULL),
  58 fErrXH(NULL),
  59 fMultHa(NULL),
  60 fErrXHa(NULL),
  61 fDistH(NULL)
  62 {
  63   // Default Constructor
  64   Reset();
  65
  66   // Histograms initialization
  67   const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
  68   const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
  69   fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
  70                        2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
  71                        2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
  72   fVertexXY->SetXTitle("X , cm");
  73   fVertexXY->SetYTitle("Y , cm");
  74   fVertexXY->SetOption("colz");
  75   fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
  76                        2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
  77   fVertexZ->SetXTitle("Z , cm");
  78   fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
  79   fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
  80   fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
  81   fErrXHa->SetLineColor(kRed);
  82   fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
  83   fMultHa->SetLineColor(kRed);
  84   fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
  85   fClu0 = new UInt_t [fgkMaxNumOfEvents];
  86   for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
  87   fAccEvents.ResetAllBits(kTRUE);
  88 }
  89
  90 //______________________________________________________________________
  91 void AliITSMeanVertexer::Reset(){
  92   // Reset averages
  93   for(Int_t i=0;i<3;i++){
  94     fWeighPosSum[i] = 0.;
  95     fWeighSigSum[i] = 0.;
  96     fAverPosSum[i] = 0.;
  97     fWeighPos[i] = 0.;
  98     fWeighSig[i] = 0.;
  99     fAverPos[i] = 0.;
 100     for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
 101     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
 102     fNoEventsContr=0;
 103     fTotContributors = 0.;
 104     if(fVertexXY)fVertexXY->Reset();
 105     if(fVertexZ)fVertexZ->Reset();
 106   }
 107 }
 108 //______________________________________________________________________
 109 Bool_t AliITSMeanVertexer::Init() {
 110   // Initialize filters
 111   // Initialize geometry
 112   // Initialize ITS classes
 113
 114   AliGeomManager::LoadGeometry();
 115   if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
 116
 117   AliITSInitGeometry initgeom;
 118   AliITSgeom *geom = initgeom.CreateAliITSgeom();
 119   if (!geom) return kFALSE;
 120   printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
 121
 122   fDetTypeRec = new AliITSDetTypeRec();
 123   fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
 124   fDetTypeRec->SetITSgeom(geom);
 125   fDetTypeRec->SetDefaults();
 126   fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
 127
 128   // Initialize filter values to their defaults
 129   SetFilterOnContributors();
 130   SetCutOnErrX();
 131   SetCutOnR();
 132   SetCutOnCls();
 133
 134   // Instatiate vertexer
 135   if (!fMode) {
 136     fVertexer = new AliITSVertexer3DTapan(1000);
 137   }
 138   else {
 139     fVertexer = new AliITSVertexer3D();
 140     fVertexer->SetDetTypeRec(fDetTypeRec);
 141     AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
 142     alias->SetWideFiducialRegion(40.,0.5);
 143     alias->SetNarrowFiducialRegion(0.5,0.5);
 144     alias->SetDeltaPhiCuts(0.5,0.025);
 145     alias->SetDCACut(0.1);
 146     alias->SetPileupAlgo(3);
 147     fVertexer->SetComputeMultiplicity(kFALSE);
 148   }
 149   return kTRUE;
 150 }
 151
 152 //______________________________________________________________________
 153 AliITSMeanVertexer::~AliITSMeanVertexer() {
 154   // Destructor
 155   delete fDetTypeRec;
 156   delete fVertexXY;
 157   delete fVertexZ;
 158   delete fMultH;
 159   delete fErrXH;
 160   delete fMultHa;
 161   delete fErrXHa;
 162   delete fDistH;
 163   delete fVertexer;
 164 }
 165
 166 //______________________________________________________________________
 167 Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
 168   // Performs SPD local reconstruction
 169   // and vertex finding
 170   // returns true in case a vertex is found
 171
 172   // Run SPD cluster finder
 173   static Int_t evcount = -1;
 174   if(evcount <0){
 175     evcount++;
 176     AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
 177     AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
 178   }
 179 //  AliCodeTimerAuto("",0);
 180   AliITSRecPointContainer::Instance()->PrepareToRead();
 181   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
 182   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
 183
 184   Bool_t vtxOK = kFALSE;
 185   UInt_t nl1=0;
 186   AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
 187   if (!fMode) {
 188     if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
 189   }
 190   else {
 191     AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
 192     nl1=rpcont->GetNClustersInLayerFast(1);
 193     if(Filter(vtx,nl1)) vtxOK = kTRUE;
 194     /*    if(vtx){
 195       if(vtxOK){
 196         printf("The vertex is OK\n");
 197       }
 198       else {
 199         printf("The vertex is NOT OK\n");
 200       }
 201       vtx->PrintStatus();
 202     }
 203     else {
 204       printf("The vertex was not reconstructed\n");
 205       }  */
 206   }
 207   delete clustersTree;
 208   if (vtxOK && fMode){
 209     new(fVertArray[fIndex])AliESDVertex(*vtx);
 210     fClu0[fIndex]=nl1;
 211     fAccEvents.SetBitNumber(fIndex);
 212     fIndex++;
 213     if(fIndex>=fgkMaxNumOfEvents){
 214       if(ComputeMean(kFALSE)){
 215         if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
 216       }
 217     }
 218   }
 219   if (vtx) delete vtx;
 220
 221   return vtxOK;
 222 }
 223
 224 //______________________________________________________________________
 225 void AliITSMeanVertexer::WriteVertices(const char *filename){
 226   // Compute mean vertex and
 227   // store it along with the histograms
 228   // in a file
 229
 230   TFile fmv(filename,"update");
 231   Bool_t itisOK = kFALSE;
 232   if(ComputeMean(kFALSE)){
 233     if(ComputeMean(kTRUE)){
 234       Double_t cov[6];
 235       cov[0] =  fAverPosSq[0][0];  // variance x
 236       cov[1] =  fAverPosSq[0][1];  // cov xy
 237       cov[2] =  fAverPosSq[1][1];  // variance y
 238       cov[3] =  fAverPosSq[0][2];  // cov xz
 239       cov[4] =  fAverPosSq[1][2];  // cov yz
 240       cov[5] =  fAverPosSq[2][2];  // variance z
 241       AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
 242       mv.SetTitle("Mean Vertex");
 243       mv.SetName("MeanVertex");
 244       // we have to add chi2 here
 245       AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
 246
 247       mv.Write(mv.GetName(),TObject::kOverwrite);
 248       vtx.Write(vtx.GetName(),TObject::kOverwrite);
 249       itisOK = kTRUE;
 250     }
 251   }
 252
 253   if(!itisOK){
 254     AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
 255   }
 256
 257   fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
 258   fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
 259   fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
 260   fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
 261   fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
 262   fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
 263   fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
 264   fmv.Close();
 265 }
 266
 267 //______________________________________________________________________
 268 Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
 269   // Apply selection criteria to events
 270   Bool_t status = kFALSE;
 271   if(!vert)return status;
 272   // Remove vertices reconstructed with vertexerZ
 273   if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
 274   Int_t ncontr = vert->GetNContributors();
 275   AliDebug(1,Form("Number of contributors = %d",ncontr));
 276   Double_t ex = vert->GetXRes();
 277   fMultH->Fill(mult);
 278   fErrXH->Fill(ex*1000.);
 279   if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
 280      ((ex*1000.)<fErrXCut)) {
 281     status = kTRUE;
 282     fMultHa->Fill(mult);
 283     fErrXHa->Fill(ex*1000.);
 284   }
 285   return status;
 286 }
 287
 288 //______________________________________________________________________
 289 void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
 290   // update mean vertex
 291   Double_t currentPos[3],currentSigma[3];
 292   vert->GetXYZ(currentPos);
 293   vert->GetSigmaXYZ(currentSigma);
 294   Bool_t goon = kTRUE;
 295   for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
 296   if(!goon)return;
 297   for(Int_t i=0;i<3;i++){
 298     fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
 299     fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
 300     fAverPosSum[i]+=currentPos[i];
 301   }
 302   for(Int_t i=0;i<3;i++){
 303     for(Int_t j=i;j<3;j++){
 304       fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
 305     }
 306   }
 307
 308   fTotContributors+=vert->GetNContributors();
 309   fVertexXY->Fill(currentPos[0],currentPos[1]);
 310   fVertexZ->Fill(currentPos[2]);
 311
 312   fNoEventsContr++;
 313 }
 314
 315 //______________________________________________________________________
 316  Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
 317    // compute mean vertex
 318    // called once with killOutliers = kFALSE and then again with
 319    // killOutliers = kTRUE
 320
 321    static Bool_t firstime = kTRUE;
 322    if(fNoEventsContr>2 && !killOutliers)return kTRUE;  //ComputeMean is
 323                              // executed only once with argument kFALSE
 324    Double_t wpos[3];
 325    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
 326
 327    if(killOutliers && firstime){
 328      firstime = kFALSE;
 329      Reset();
 330    }
 331
 332    for(Int_t i=0;i<fVertArray.GetEntries();i++){
 333      AliESDVertex* vert = NULL;
 334
 335      // second pass
 336      if(killOutliers && fAccEvents.TestBitNumber(i)){
 337        vert=(AliESDVertex*)fVertArray[i];
 338        Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
 339        dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
 340        dist=sqrt(dist)*10.;    // distance in mm
 341        fDistH->Fill(dist);
 342        if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
 343      }
 344
 345      if(!fAccEvents.TestBitNumber(i))continue;
 346      vert=(AliESDVertex*)fVertArray[i];
 347      AddToMean(vert);
 348    }
 349    if(fNoEventsContr < 2) return kFALSE;
 350    Double_t nevents = fNoEventsContr;
 351    for(Int_t i=0;i<3;i++){
 352      fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
 353      fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
 354      fAverPos[i] = fAverPosSum[i]/nevents;
 355    }
 356    for(Int_t i=0;i<3;i++){
 357      for(Int_t j=i;j<3;j++){
 358        fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
 359        fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
 360      }
 361    }
 362    if(killOutliers)ResetArray();
 363    fAverContributors = fTotContributors/nevents;
 364    return kTRUE;
 365  }