ITS/AliITSVertex.cxx

   1 #include "stdlib.h"
   2 #include <TMath.h>
   3 #include <TRandom.h>
   4 #include <TObjArray.h>
   5 #include <TROOT.h>
   6 #include <TFile.h>
   7 #include <TTree.h>
   8 #include <iostream.h>
   9
  10 #include "AliRun.h"
  11 #include "AliITS.h"
  12 #include "AliITSgeom.h"
  13 #include "AliITSRecPoint.h"
  14 #include "AliGenerator.h"
  15 #include "AliMagF.h"
  16
  17 #include <TH1.h>
  18 #include <TF1.h>
  19 #include <TCanvas.h>
  20 #include <AliITSVertex.h>
  21 #include <TObjArray.h>
  22 #include <TObject.h>
  23
  24
  25 ClassImp(AliITSVertex)
  26
  27 //////////////////////////////////////////////////////////////////////
  28 // AliITSVertex is a class for full 3D primary vertex finding       //
  29 //                                                                  //                                                           //
  30 // Version: 1                                                       //                                   //
  31 // Written by Giuseppe Lo Re and Francesco Riggi                    //
  32 // Giuseppe.Lore@ct.infn.it                                         //                           //
  33 // Franco.Riggi@ct.infn.it                                          //                           //
  34 // Release date: May 2001                                           //                                                                                           //
  35 //                                                                  //                                                   //
  36 //                                                                  //                                                   //
  37 //////////////////////////////////////////////////////////////////////
  38
  39
  40 //______________________________________________________________________________
  41
  42 AliITSVertex::AliITSVertex() {
  43
  44    fPosition = new Double_t[3];
  45    fResolution = new Double_t[3];
  46    fSNR = new Double_t[3];
  47
  48    AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS");
  49    AliITSgeom *g2 = ((AliITS*)aliits)->GetITSgeom();
  50
  51    TClonesArray  *recpoints = aliits->RecPoints();
  52    AliITSRecPoint *pnt;
  53
  54    Int_t NoPoints1 = 0;
  55    Int_t NoPoints2 = 0;
  56    Double_t Vzero[3];
  57    Double_t AsPar[2];
  58
  59 //------------ Rough Vertex evaluation ---------------------------------
  60
  61    Int_t i,npoints,ipoint,j,k,max,BinMax;
  62    Double_t ZCentroid;
  63    Float_t l[3], p[3];
  64
  65    Double_t mxpiu = 0;
  66    Double_t mxmeno = 0;
  67    Double_t mypiu = 0;
  68    Double_t mymeno = 0;
  69
  70    TH1F *hITSz1         = new TH1F("hITSz1","",100,-14.35,14.35);
  71
  72    for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
  73    {
  74            aliits->ResetRecPoints();
  75        gAlice->TreeR()->GetEvent(i);
  76            npoints = recpoints->GetEntries();
  77
  78            for (ipoint=0;ipoint<npoints;ipoint++) {
  79                pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
  80                l[0]=pnt->GetX();
  81                l[1]=0;
  82                l[2]=pnt->GetZ();
  83                    g2->LtoG(i, l, p);
  84                    if(i<80 && TMath::Abs(p[2])<14.35) {
  85                         if(p[0]>0) mxpiu++;
  86                         if(p[0]<0) mxmeno++;
  87                         if(p[1]>0) mypiu++;
  88                         if(p[1]<0) mymeno++;
  89                    hITSz1->Fill(p[2]);
  90                }
  91                    if(i>=80 && TMath::Abs(p[2])<14.35) NoPoints2++;
  92        }
  93    }
  94
  95    NoPoints1 = (Int_t)(hITSz1->GetEntries());
  96
  97    AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
  98    AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno);
  99
 100    Vzero[0] = 5.24441*AsPar[0];
 101    Vzero[1] = 5.24441*AsPar[1];
 102
 103    ZCentroid=  TMath::Abs(hITSz1->GetMean());
 104    Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
 105                   -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)
 106               -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::Power(ZCentroid,6);
 107
 108    if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
 109
 110    cout << "\nXvzero: " << Vzero[0] << " cm" << "";
 111    cout << "\nYvzero: " << Vzero[1] << " cm" << "";
 112    cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";
 113
 114    delete hITSz1;
 115
 116    Double_t dphi,r,DeltaZ,a,b;
 117    Int_t np1=0;
 118    Int_t np2=0;
 119    Int_t niter=0;
 120
 121    Double_t DeltaPhiZ = 0.08;
 122
 123    Float_t B[3];
 124    Float_t origin[3];
 125    for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
 126    gAlice->Field()->Field(origin,B);
 127
 128    DeltaPhiZ = DeltaPhiZ*B[2]/2;
 129    Double_t DeltaPhiXY = 1.0;
 130
 131 //   cout << "\nDeltaPhiZ: " << DeltaPhiZ << " deg" << "\n";
 132 //   cout << "DeltaPhiXY: " << DeltaPhiXY << " deg" << "\n";
 133
 134    Double_t *Z1, *Z2, *Y1, *Y2, *X1, *X2, *phi1, *phi2, *r1, *r2;
 135    Z1=new Double_t[NoPoints1];
 136    Z2=new Double_t[NoPoints2];
 137    Y1=new Double_t[NoPoints1];
 138    Y2=new Double_t[NoPoints2];
 139    X1=new Double_t[NoPoints1];
 140    X2=new Double_t[NoPoints2];
 141    phi1=new Double_t[NoPoints1];
 142    phi2=new Double_t[NoPoints2];
 143    r1=new Double_t[NoPoints1];
 144    r2=new Double_t[NoPoints2];
 145
 146    DeltaZ = 4.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2);
 147    Float_t MulFactorZ = 28000./(Float_t)NoPoints1;
 148    Int_t nbin=(Int_t)((DeltaZ/0.005)/MulFactorZ);
 149    Int_t nbinxy=250;
 150    Int_t *VectorBinZ,*VectorBinXY;
 151    VectorBinZ=new Int_t[nbin];
 152    VectorBinXY=new Int_t[nbinxy];
 153    Float_t f1= 0;
 154    Float_t f2= 0;
 155    Double_t sigma,MediaFondo;
 156
 157    TH1D *hITSZv         = new TH1D("hITSZv","",nbin,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);
 158    TH1D *hITSXv         = new TH1D("hITSXv","",nbinxy,-3,3);
 159    TH1D *hITSYv         = new TH1D("hITSYv","",nbinxy,-3,3);
 160
 161 //   cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";
 162
 163
 164    start:
 165
 166    hITSZv->Add(hITSZv,-1.);
 167    hITSXv->Add(hITSXv,-1.);
 168    hITSYv->Add(hITSYv,-1.);
 169
 170    np1=np2=0;
 171
 172
 173
 174    for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
 175    {
 176         aliits->ResetRecPoints();
 177         gAlice->TreeR()->GetEvent(i);
 178            npoints = recpoints->GetEntries();
 179            for (ipoint=0;ipoint<npoints;ipoint++) {
 180
 181                     pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
 182                 l[0]=pnt->GetX();
 183                 l[1]=0;
 184                 l[2]=pnt->GetZ();
 185                 g2->LtoG(i, l, p);
 186
 187                         if(i<80 && TMath::Abs(p[2])<14.35) {
 188                     p[0]=p[0]-Vzero[0];
 189                     p[1]=p[1]-Vzero[1];
 190                 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 191                         Y1[np1]=p[1];
 192                         X1[np1]=p[0];
 193                         Z1[np1]=p[2];
 194                         r1[np1]=r;
 195                         phi1[np1]=PhiFunc(p);
 196                         np1++;
 197                         }
 198
 199                         if(i>=80 &&  TMath::Abs(p[2])<14.35) {
 200                     p[0]=p[0]-Vzero[0];
 201                     p[1]=p[1]-Vzero[1];
 202                 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 203                         Y2[np2]=p[1];
 204                         X2[np2]=p[0];
 205                         Z2[np2]=p[2];
 206                         r2[np2]=r;
 207                         phi2[np2]=PhiFunc(p);
 208                         np2++;
 209                         }
 210
 211           }
 212    }
 213
 214 //------------------ Correlation between rec points ----------------------
 215
 216    cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl;
 217
 218    for(j=0; j<(np2)-1; j++) {
 219          for(k=0; k<(np1)-1; k++) {
 220                 dphi=TMath::Abs(phi1[k]-phi2[j]);
 221                 if(dphi>180) dphi = 360-dphi;
 222                 if(dphi<DeltaPhiZ && TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Vzero[2])
 223                          <DeltaZ) hITSZv->Fill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1));
 224          }
 225    }
 226
 227    //cout << "\nNumber of used pairs: \n";
 228    //cout << hITSZv->GetEntries() << '\n' << '\n';
 229    a = Vzero[2]-DeltaZ;
 230    b = Vzero[2]+DeltaZ;
 231    max=(Int_t) hITSZv->GetMaximum();
 232    BinMax=hITSZv->GetMaximumBin();
 233    sigma=0;
 234    for(i=0;i<nbin;i++) VectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
 235    for(i=0;i<10;i++) f1=f1+VectorBinZ[i]/10;
 236    for(i=nbin-10;i<nbin;i++) f2=f2+VectorBinZ[i]/10;
 237    MediaFondo=(f1+f2)/2;
 238    for(i=0;i<nbin;i++) {
 239           if(VectorBinZ[i]-MediaFondo>(max-MediaFondo)*0.4 &&
 240           VectorBinZ[i]-MediaFondo<(max-MediaFondo)*0.7) {
 241                  sigma=hITSZv->GetBinCenter(BinMax)-hITSZv->GetBinCenter(i);
 242                  sigma=TMath::Abs(sigma);
 243                  if(sigma==0) sigma=0.05;
 244       }
 245    }
 246
 247    cout << "f1 " <<f1 <<endl;
 248    cout << "f2 " <<f2 <<endl;
 249    cout << "GetMaximumBin " <<hITSZv->GetMaximumBin() <<endl;
 250    cout << "nbin " <<nbin <<endl;
 251    cout << "max " << hITSZv->GetBinContent(BinMax)<<endl;
 252    cout << "sigma " <<sigma <<endl;
 253    cout << "Fondo " <<   MediaFondo<< endl;
 254
 255    TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);
 256
 257    fz->SetParameter(0,max);
 258    if(niter==0) {Double_t temp = hITSZv->GetBinCenter(BinMax); Vzero[2]=temp;}
 259    fz->SetParameter(1,Vzero[2]);
 260    fz->SetParameter(2,sigma);
 261    fz->SetParameter(3,MediaFondo);
 262    fz->SetParLimits(2,0,999);
 263    fz->SetParLimits(3,0,999);
 264
 265    hITSZv->Fit("fz","RMEQ0");
 266
 267    fSNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
 268    if(fSNR[2]<0.) {
 269      cout << "\nNegative Signal to noise ratio for z!!!" << endl;
 270      cout << "The algorithm cannot find the z vertex position." << endl;
 271      exit(123456789);
 272    }
 273    else
 274    {
 275      fPosition[2] = fz->GetParameter(1);
 276      if(fPosition[2]<0)
 277        {
 278          fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000;
 279        }
 280      else
 281        {
 282          fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000;
 283        }
 284    }
 285    fResolution[2] = fz->GetParError(1);
 286
 287
 288
 289    for(j=0; j<(np2)-1; j++) {
 290          for(k=0; k<(np1)-1; k++) {
 291                 dphi=TMath::Abs(phi1[k]-phi2[j]);
 292                 if(dphi>180) dphi = 360-dphi;
 293                 if(dphi>DeltaPhiXY) continue;
 294             if(TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-fPosition[2])
 295                 <4*fResolution[2]) {
 296                    hITSXv->Fill(Vzero[0]+(X2[j]-((X2[j]-X1[k])/(Y2[j]-Y1[k]))*Y2[j]));
 297                    hITSYv->Fill(Vzero[1]+(Y2[j]-((Y2[j]-Y1[k])/(X2[j]-X1[k]))*X2[j]));
 298             }
 299          }
 300    }
 301
 302    TF1 *fx = new TF1
 303    ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[0]-0.5,Vzero[0]+0.5);
 304
 305    max=(Int_t) hITSXv->GetMaximum();
 306    BinMax=hITSXv->GetMaximumBin();
 307    sigma=0;
 308    f1=f2=0;
 309    for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
 310    for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
 311    for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
 312    MediaFondo=(f1+f2)/2;
 313    for(i=0;i<nbinxy;i++) {
 314           if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 &&
 315           VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
 316                  sigma=hITSXv->GetBinCenter(BinMax)-hITSXv->GetBinCenter(i);
 317                  sigma=TMath::Abs(sigma);
 318                  if(sigma==0) sigma=0.05;
 319       }
 320    }
 321
 322    /*cout << "f1 " <<f1 <<endl;
 323    cout << "f2 " <<f2 <<endl;
 324    cout << "GetMaximumBin " <<hITSXv->GetMaximumBin() <<endl;
 325    cout << "max " << hITSXv->GetBinContent(BinMax)<<endl;
 326    cout << "sigma " <<sigma <<endl;
 327    cout << "Fondo " <<   MediaFondo<< endl;
 328    cout << "nbinxy " <<nbinxy <<endl;*/
 329
 330    fx->SetParameter(0,max);
 331    fx->SetParameter(1,Vzero[0]);
 332    fx->SetParameter(2,sigma);
 333    fx->SetParameter(3,MediaFondo);
 334
 335    hITSXv->Fit("fx","RMEQ0");
 336
 337    fSNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
 338    if(fSNR[0]<0.) {
 339      cout << "\nNegative Signal to noise ratio for x!!!" << endl;
 340      cout << "The algorithm cannot find the x vertex position." << endl;
 341      exit(123456789);
 342    }
 343    else
 344    {
 345    fPosition[0]=fx->GetParameter(1);
 346    fResolution[0]=fx->GetParError(1);
 347    }
 348
 349    TF1 *fy = new TF1
 350    ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[1]-0.5,Vzero[1]+0.5);
 351
 352    max=(Int_t) hITSYv->GetMaximum();
 353    BinMax=hITSYv->GetMaximumBin();
 354    sigma=0;
 355    f1=f2=0;
 356    for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
 357    for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
 358    for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
 359    MediaFondo=(f1+f2)/2;
 360    for(i=0;i<nbinxy;i++) {
 361           if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 &&
 362           VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
 363                  sigma=hITSYv->GetBinCenter(BinMax)-hITSYv->GetBinCenter(i);
 364                  sigma=TMath::Abs(sigma);
 365                  if(sigma==0) sigma=0.05;
 366       }
 367    }
 368
 369    /*cout << "f1 " <<f1 <<endl;
 370    cout << "f2 " <<f2 <<endl;
 371    cout << "GetMaximumBin " <<hITSYv->GetMaximumBin() <<endl;
 372    cout << "max " << hITSYv->GetBinContent(BinMax)<<endl;
 373    cout << "sigma " <<sigma <<endl;
 374    cout << "Fondo " <<   MediaFondo<< endl;
 375    cout << "nbinxy " <<nbinxy <<endl;*/
 376
 377    fy->SetParameter(0,max);
 378    fy->SetParameter(1,Vzero[1]);
 379    fy->SetParameter(2,sigma);
 380    fy->SetParameter(3,MediaFondo);
 381
 382    hITSYv->Fit("fy","RMEQ0");
 383
 384    fSNR[1] = fy->GetParameter(0)/fy->GetParameter(3);
 385    if(fSNR[1]<0.) {
 386      cout << "\nNegative Signal to noise ratio for y!!!" << endl;
 387      cout << "The algorithm cannot find the y vertex position." << endl;
 388      exit(123456789);
 389    }
 390    else
 391    {
 392    fPosition[1]=fy->GetParameter(1);
 393    fResolution[1]=fy->GetParError(1);
 394    }
 395
 396    niter++;
 397
 398    Vzero[0] = fPosition[0];
 399    Vzero[1] = fPosition[1];
 400    Vzero[2] = fPosition[2];
 401
 402    if(niter<3) goto start;  // number of iterations
 403
 404
 405    cout<<"Number of iterations: "<<niter<<endl;
 406
 407    delete [] Z1;
 408    delete [] Z2;
 409    delete [] Y1;
 410    delete [] Y2;
 411    delete [] X1;
 412    delete [] X2;
 413    delete [] r1;
 414    delete [] r2;
 415    delete [] phi1;
 416    delete [] phi2;
 417    delete [] VectorBinZ;
 418    delete [] VectorBinXY;
 419
 420    delete hITSZv;
 421    delete hITSXv;
 422    delete hITSYv;
 423
 424 }
 425
 426
 427 //______________________________________________________________________________
 428
 429 AliITSVertex::~AliITSVertex() {
 430    delete [] fPosition;
 431    delete [] fResolution;
 432    delete [] fSNR;
 433 }
 434
 435 //______________________________________________________________________________
 436
 437 Double_t AliITSVertex::PhiFunc(Float_t p[]) {
 438          Double_t phi=0;
 439          if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
 440          if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
 441          if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
 442          if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
 443          return phi;
 444 }
 445
 446 //______________________________________________________________________________
 447
 448