ITS/AliITSVertexerIons.cxx

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