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