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