ITS/AliITSVertexerZ.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 #include "AliITSVertexerZ.h"
  16 #include<TBranch.h>
  17 #include<TClonesArray.h>
  18 #include<TFile.h>
  19 #include<TH1.h>
  20 #include <TString.h>
  21 #include<TTree.h>
  22 #include "AliITSLoader.h"
  23 #include "AliITSgeom.h"
  24 #include "AliITSDetTypeRec.h"
  25 #include "AliITSRecPoint.h"
  26 #include "AliITSZPoint.h"
  27 #include "AliHeader.h"
  28 #include "AliGenEventHeader.h"
  29
  30 /////////////////////////////////////////////////////////////////
  31 // this class implements a fast method to determine
  32 // the Z coordinate of the primary vertex
  33 // for p-p collisions it seems to give comparable or better results
  34 // with respect to what obtained with AliITSVertexerPPZ
  35 // It can be used successfully with Pb-Pb collisions
  36 ////////////////////////////////////////////////////////////////
  37
  38 ClassImp(AliITSVertexerZ)
  39
  40
  41
  42 //______________________________________________________________________
  43 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
  44 fFirstL1(0),
  45 fLastL1(0),
  46 fFirstL2(0),
  47 fLastL2(0),
  48 fDiffPhiMax(0),
  49 fX0(0.),
  50 fY0(0.),
  51 fZFound(0),
  52 fZsig(0.),
  53 fZCombc(0),
  54 fLowLim(0.),
  55 fHighLim(0.),
  56 fStepCoarse(0),
  57 fTolerance(0.),
  58 fMaxIter(0),
  59 fWindowWidth(0) {
  60   // Default constructor
  61   SetDiffPhiMax();
  62   SetFirstLayerModules();
  63   SetSecondLayerModules();
  64   SetLowLimit();
  65   SetHighLimit();
  66   SetBinWidthCoarse();
  67   SetTolerance();
  68   SetPPsetting();
  69   ConfigIterations();
  70   SetWindowWidth();
  71 }
  72
  73 //______________________________________________________________________
  74 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn),
  75 fFirstL1(0),
  76 fLastL1(0),
  77 fFirstL2(0),
  78 fLastL2(0),
  79 fDiffPhiMax(0),
  80 fX0(x0),
  81 fY0(y0),
  82 fZFound(0),
  83 fZsig(0.),
  84 fZCombc(0),
  85 fLowLim(0.),
  86 fHighLim(0.),
  87 fStepCoarse(0),
  88 fTolerance(0.),
  89 fMaxIter(0),
  90 fWindowWidth(0) {
  91   // Standard Constructor
  92   SetDiffPhiMax();
  93   SetFirstLayerModules();
  94   SetSecondLayerModules();
  95   SetLowLimit();
  96   SetHighLimit();
  97   SetBinWidthCoarse();
  98   SetTolerance();
  99   SetPPsetting();
 100   ConfigIterations();
 101   SetWindowWidth();
 102
 103 }
 104
 105 //______________________________________________________________________
 106 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr),
 107 fFirstL1(vtxr.fFirstL1),
 108 fLastL1(vtxr.fLastL1),
 109 fFirstL2(vtxr.fFirstL2),
 110 fLastL2(vtxr.fLastL2),
 111 fDiffPhiMax(vtxr.fDiffPhiMax),
 112 fX0(vtxr.fX0),
 113 fY0(vtxr.fY0),
 114 fZFound(vtxr.fZFound),
 115 fZsig(vtxr.fZsig),
 116 fZCombc(vtxr.fZCombc),
 117 fLowLim(vtxr.fLowLim),
 118 fHighLim(vtxr.fHighLim),
 119 fStepCoarse(vtxr.fStepCoarse),
 120 fTolerance(vtxr.fTolerance),
 121 fMaxIter(vtxr.fMaxIter),
 122 fWindowWidth(vtxr.fWindowWidth){
 123   // Copy constructor
 124
 125 }
 126
 127 //______________________________________________________________________
 128 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ&  vtxr ){
 129   // Assignment operator
 130
 131   this->~AliITSVertexerZ();
 132   new(this) AliITSVertexerZ(vtxr);
 133   return *this;
 134 }
 135
 136 //______________________________________________________________________
 137 AliITSVertexerZ::~AliITSVertexerZ() {
 138   // Destructor
 139   delete fZCombc;
 140 }
 141
 142 //______________________________________________________________________
 143 void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
 144   // configure the iterative procedure to gain efficiency for
 145   // pp events with very low multiplicity
 146   Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
 147   fMaxIter=noiter;
 148   if(noiter>5){
 149     Error("ConfigIterations","Maximum number of iterations is 5\n");
 150     fMaxIter=5;
 151   }
 152   for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
 153   if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
 154 }
 155
 156 //______________________________________________________________________
 157 Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax) const {
 158   // Finds a region around a peak in the Z histogram
 159   // Case of 2 peaks is treated
 160   Int_t imax=h->GetNbinsX();
 161   Float_t maxval=0;
 162   Int_t bi1=h->GetMaximumBin();
 163   Int_t bi2=0;
 164   for(Int_t i=imax;i>=1;i--){
 165     if(h->GetBinContent(i)>maxval){
 166       maxval=h->GetBinContent(i);
 167       bi2=i;
 168     }
 169   }
 170   Int_t npeaks=0;
 171
 172   if(bi1==bi2){
 173     binmin=bi1-3;
 174     binmax=bi1+3;
 175     npeaks=1;
 176   }else{
 177     TH1F *copy = new TH1F(*h);
 178     copy->SetBinContent(bi1,0.);
 179     copy->SetBinContent(bi2,0.);
 180     Int_t l1=TMath::Max(bi1-3,1);
 181     Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
 182     Float_t cont1=copy->Integral(l1,l2);
 183     Int_t ll1=TMath::Max(bi2-3,1);
 184     Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
 185     Float_t cont2=copy->Integral(ll1,ll2);
 186     if(cont1>cont2){
 187       binmin=l1;
 188       binmax=l2;
 189       npeaks=1;
 190     }
 191     if(cont2>cont1){
 192       binmin=ll1;
 193       binmax=ll2;
 194       npeaks=1;
 195     }
 196     if(cont1==cont2){
 197       binmin=l1;
 198       binmax=ll2;
 199       if(bi2-bi1==1) npeaks=1;
 200       else npeaks=2;
 201     }
 202     delete copy;
 203   }
 204   return npeaks;
 205 }
 206 //______________________________________________________________________
 207 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
 208   // Defines the AliESDVertex for the current event
 209   VertexZFinder(evnumber);
 210   Int_t ntrackl=0;
 211   for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
 212     if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
 213     if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){
 214       Float_t diffPhiMaxOrig=fDiffPhiMax;
 215       fDiffPhiMax=GetPhiMaxIter(iteraz);
 216       VertexZFinder(evnumber);
 217       fDiffPhiMax=diffPhiMaxOrig;
 218     }
 219   }
 220   FindMultiplicity(evnumber);
 221   return fCurrentVertex;
 222 }
 223
 224
 225
 226
 227 //______________________________________________________________________
 228 void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
 229   // Defines the AliESDVertex for the current event
 230   fCurrentVertex = 0;
 231   AliRunLoader *rl =AliRunLoader::GetRunLoader();
 232   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
 233   AliITSgeom* geom = itsLoader->GetITSgeom();
 234   itsLoader->LoadRecPoints();
 235   rl->GetEvent(evnumber);
 236
 237   AliITSDetTypeRec detTypeRec;
 238
 239   TTree *tR = itsLoader->TreeR();
 240   detTypeRec.SetTreeAddressR(tR);
 241   TClonesArray *itsRec  = 0;
 242   // lc and gc are local and global coordinates for layer 1
 243   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
 244   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
 245   // lc2 and gc2 are local and global coordinates for layer 2
 246   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
 247   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
 248
 249   itsRec = detTypeRec.RecPoints();
 250   TBranch *branch;
 251   branch = tR->GetBranch("ITSRecPoints");
 252
 253   Int_t nrpL1 = 0;
 254   Int_t nrpL2 = 0;
 255   // By default fFirstL1=0 and fLastL1=79
 256   // This loop counts the number of recpoints on layer1 (central modules)
 257   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 258     // Keep only central modules
 259     //    if(module%4==0 || module%4==3)continue;
 260     //   cout<<"Procesing module "<<module<<" ";
 261     branch->GetEvent(module);
 262     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
 263     nrpL1+= itsRec->GetEntries();
 264     detTypeRec.ResetRecPoints();
 265   }
 266   //By default fFirstL2=80 and fLastL2=239
 267   //This loop counts the number of RP on layer 2
 268   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 269     branch->GetEvent(module);
 270     nrpL2+= itsRec->GetEntries();
 271     detTypeRec.ResetRecPoints();
 272   }
 273   if(nrpL1 == 0 || nrpL2 == 0){
 274     ResetHistograms();
 275     itsLoader->UnloadRecPoints();
 276     fCurrentVertex = new AliESDVertex(0.,5.3,-2);
 277     return;
 278   }
 279   // The vertex finding is attempted only if the number of RP is !=0 on
 280   // both layers
 281   Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
 282   Float_t *yc1 = new Float_t [nrpL1];
 283   Float_t *zc1 = new Float_t [nrpL1];
 284   Float_t *phi1 = new Float_t [nrpL1];
 285   Float_t *err1 = new Float_t [nrpL1];
 286   Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
 287   Float_t *yc2 = new Float_t [nrpL2];
 288   Float_t *zc2 = new Float_t [nrpL2];
 289   Float_t *phi2 = new Float_t [nrpL2];
 290   Float_t *err2 = new Float_t [nrpL2];
 291   Int_t ind = 0;// running index for RP
 292   // Force a coarse bin size of 200 microns if the number of clusters on layer 2
 293   // is low
 294   if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
 295   // By default nbincoarse=(10+10)/0.01=2000
 296   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
 297   if(fZCombc)delete fZCombc;
 298   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
 299
 300   // Loop on modules of layer 1
 301
 302   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 303     //    if(module%4==0 || module%4==3)continue;
 304     branch->GetEvent(module);
 305     Int_t nrecp1 = itsRec->GetEntries();
 306     for(Int_t j=0;j<nrecp1;j++){
 307       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 308       // Local coordinates of this recpoint
 309       lc[0]=recp->GetDetLocalX();
 310       lc[2]=recp->GetDetLocalZ();
 311       geom->LtoG(module,lc,gc);
 312       // Global coordinates of this recpoints
 313       gc[0]-=fX0; // Possible beam offset in the bending plane
 314       gc[1]-=fY0; //   "               "
 315       xc1[ind]=gc[0];
 316       yc1[ind]=gc[1];
 317       zc1[ind]=gc[2];
 318       // azimuthal angle is computed in the interval 0 --> 2*pi
 319       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
 320       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
 321       err1[ind]=recp->GetSigmaZ2();
 322       ind++;
 323     }
 324     detTypeRec.ResetRecPoints();
 325   }
 326   ind = 0; // the running index is reset for Layer 2
 327   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 328     branch->GetEvent(module);
 329     Int_t nrecp2 = itsRec->GetEntries();
 330     for(Int_t j=0;j<nrecp2;j++){
 331       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 332       lc[0]=recp->GetDetLocalX();
 333       lc[2]=recp->GetDetLocalZ();
 334       geom->LtoG(module,lc,gc);
 335       gc[0]-=fX0;
 336       gc[1]-=fY0;
 337       xc2[ind]=gc[0];
 338       yc2[ind]=gc[1];
 339       zc2[ind]=gc[2];
 340       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
 341       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
 342       err2[ind]=recp->GetSigmaZ2();
 343       ind++;
 344     }
 345     detTypeRec.ResetRecPoints();
 346   }
 347
 348 /* Test the ffect of mutiple scatternig on error. Negligible
 349   // Multiple scattering
 350   Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
 351   Float_t beta2=beta*beta;
 352   Float_t p2=pmed*pmed;
 353   Float_t rBP=3; //Beam Pipe radius = 3cm
 354   Float_t dBP=0.08/35.3; // 800 um of Be
 355   Float_t dL1=0.01; //approx. 1% of radiation length
 356   Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
 357   Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
 358 */
 359   TClonesArray *points = new TClonesArray("AliITSZPoint",nrpL1*nrpL2);
 360   TClonesArray &pts = *points;
 361   Int_t nopoints =0;
 362   for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
 363     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
 364     for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
 365       Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
 366       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
 367       if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
 368         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
 369 //      Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope
 370         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
 371         Float_t ezr0q=(r2*r2*err1[i]+r1*r1*err2[j])/(r2-r1)/(r2-r1); //error on Z @ null radius
 372         /*
 373         ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1
 374         ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe
 375         */
 376         new(pts[nopoints++])AliITSZPoint(zr0,ezr0q);
 377
 378         fZCombc->Fill(zr0);
 379       }
 380     }
 381   }
 382   delete [] xc1;
 383   delete [] yc1;
 384   delete [] zc1;
 385   delete [] phi1;
 386   delete [] err1;
 387   delete [] xc2;
 388   delete [] yc2;
 389   delete [] zc2;
 390   delete [] phi2;
 391   delete [] err2;
 392
 393   points->Sort();
 394
 395   Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
 396   if(contents<1.){
 397     //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
 398     ResetHistograms();
 399     itsLoader->UnloadRecPoints();
 400     fCurrentVertex = new AliESDVertex(0.,5.3,-1);
 401     return;
 402   }
 403
 404   TH1F *hc = fZCombc;
 405
 406
 407   if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(3);
 408   Int_t binmin,binmax;
 409   Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);
 410   if(nPeaks==2)AliWarning("2 peaks found");
 411   Float_t zm =0.;
 412   Float_t ezm =0.;
 413   Float_t lim1 = hc->GetBinLowEdge(binmin);
 414   Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
 415
 416   if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
 417     Float_t c=(lim1+lim2)/2.;
 418     lim1=c-fWindowWidth/2.;
 419     lim2=c+fWindowWidth/2.;
 420   }
 421   Int_t niter = 0, ncontr=0;
 422   do {
 423     // symmetrization
 424     if(zm  !=0.){
 425       Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
 426       lim1=zm - semilarg;
 427       lim2=zm + semilarg;
 428     }
 429
 430     zm=0.;
 431     ezm=0.;
 432     ncontr=0;
 433     for(Int_t i =0; i<points->GetEntriesFast(); i++){
 434       AliITSZPoint* p=(AliITSZPoint*)points->UncheckedAt(i);
 435       if(p->GetZ()>lim1 && p->GetZ()<lim2){
 436         Float_t deno = p->GetErrZ();
 437         zm+=p->GetZ()/deno;
 438         ezm+=1./deno;
 439         ncontr++;
 440       }
 441     }
 442     zm/=ezm;
 443     ezm=TMath::Sqrt(1./ezm);
 444     niter++;
 445   } while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
 446   fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
 447   fCurrentVertex->SetTitle("vertexer: B");
 448   ResetHistograms();
 449   itsLoader->UnloadRecPoints();
 450   return;
 451 }
 452
 453 //_____________________________________________________________________
 454 void AliITSVertexerZ::ResetHistograms(){
 455   // delete TH1 data members
 456   if(fZCombc)delete fZCombc;
 457   fZCombc = 0;
 458 }
 459
 460 //______________________________________________________________________
 461 void AliITSVertexerZ::FindVertices(){
 462   // computes the vertices of the events in the range FirstEvent - LastEvent
 463   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 464   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 465   itsLoader->ReloadRecPoints();
 466   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 467     //  cout<<"Processing event "<<i<<endl;
 468     rl->GetEvent(i);
 469     FindVertexForCurrentEvent(i);
 470     if(fCurrentVertex){
 471       WriteCurrentVertex();
 472     }
 473   }
 474 }
 475
 476 //________________________________________________________
 477 void AliITSVertexerZ::PrintStatus() const {
 478   // Print current status
 479   cout <<"=======================================================\n";
 480   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
 481   cout <<fLastL1<<endl;
 482   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
 483   cout <<fLastL2<<endl;
 484   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
 485   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
 486   cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl;
 487   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
 488   cout <<" Debug flag: "<<fDebug<<endl;
 489   cout <<"First event to be processed "<<fFirstEvent;
 490   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
 491   if(fZCombc){
 492     cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
 493   }
 494   else{
 495     cout<<"fZCombc does not exist\n";
 496   }
 497
 498   cout <<"=======================================================\n";
 499 }
 500