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