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