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