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