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