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