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 primary vertex finding.
  29 //
  30 // Version: 0
  31 // Written by Giuseppe Lo Re and Francesco Riggi
  32 // Giuseppe.Lore@ct.infn.it
  33 // Franco.Riggi@ct.infn.it
  34 // Marh 2001
  35 //
  36 //
  37 ///////////////////////////////////////////////////////////////////////
  38
  39
  40 //______________________________________________________________________________
  41
  42 AliITSVertex::AliITSVertex() {
  43
  44    fPosition = new Double_t[3];
  45    fResolution = new Double_t[3];
  46    fSNR = new Double_t[3];
  47
  48    AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS");
  49    AliITSgeom *g2 = ((AliITS*)aliits)->GetITSgeom();
  50
  51    TClonesArray  *recpoints = aliits->RecPoints();
  52    AliITSRecPoint *pnt;
  53    TRandom rnd;
  54
  55    Int_t NoPoints1 = 0;
  56    Int_t NoPoints2 = 0;
  57    Double_t Vzero[3];
  58    Double_t AsPar[2];
  59
  60 //------------Rough Vertex------------------------------------------------------
  61
  62    Int_t i,npoints,ipoint,j,k;
  63    Double_t ZCentroid;
  64    Float_t l[3], p[3];
  65
  66    TH1F *hITSx1pos      = new TH1F("hITSx1pos","",100,0.,4.2);
  67    TH1F *hITSx1neg      = new TH1F("hITSx1neg","",100,-4.2,0.);
  68    TH1F *hITSy1pos      = new TH1F("hITSy1pos","",100,0.,4.2);
  69    TH1F *hITSy1neg      = new TH1F("hITSy1neg","",100,-4.2,0.);
  70    TH1F *hITSz1         = new TH1F("hITSz1","",100,-14.35,14.35);
  71
  72    for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
  73    {
  74         aliits->ResetRecPoints();
  75         gAlice->TreeR()->GetEvent(i);
  76
  77            npoints = recpoints->GetEntries();
  78            for (ipoint=0;ipoint<npoints;ipoint++) {
  79                pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
  80                l[0]=pnt->GetX();
  81                l[1]=0;
  82                l[2]=pnt->GetZ();
  83                    g2->LtoG(i, l, p);
  84                    if(i<80 && TMath::Abs(p[2])<14.35) {
  85                         if(p[0]>0) hITSx1pos->Fill(p[0]);
  86                         if(p[0]<0) hITSx1neg->Fill(p[0]);
  87                         if(p[1]>0) hITSy1pos->Fill(p[1]);
  88                         if(p[1]<0) hITSy1neg->Fill(p[1]);
  89                    hITSz1->Fill(p[2]);
  90                }
  91                    if(i>=80 && TMath::Abs(p[2])<14.35) (NoPoints2)++;
  92        }
  93    }
  94
  95    NoPoints1 = (Int_t)(hITSz1->GetEntries());
  96
  97    Double_t mxpiu = hITSx1pos->GetEntries();
  98    Double_t mxmeno = hITSx1neg->GetEntries();
  99    Double_t mypiu = hITSy1pos->GetEntries();
 100    Double_t mymeno = hITSy1neg->GetEntries();
 101
 102    AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
 103    AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno);
 104
 105    Vzero[0] = 5.24434*AsPar[0];
 106    Vzero[1] = 5.24434*AsPar[1];
 107
 108    ZCentroid=  TMath::Abs(hITSz1->GetMean());
 109    Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
 110                   -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)
 111               -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::Power(ZCentroid,6);
 112
 113    if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
 114
 115    //cout << "\nCentroid and RMS of hIITSz1: \n";
 116    //cout << hITSz1->GetMean() << "       "  << hITSz1->GetRMS()  <<endl;
 117    //cout << "\nAsPar[0]: " << AsPar[0];
 118    //cout << "\nAsPar[1]: " << AsPar[1];
 119    //cout << "\nAsPar[2]: " << AsPar[2];
 120    cout << "\nXvzero: " << Vzero[0] << " cm" << "";
 121    cout << "\nYvzero: " << Vzero[1] << " cm" << "";
 122    cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";
 123
 124    delete hITSx1pos;
 125    delete hITSx1neg;
 126    delete hITSy1pos;
 127    delete hITSy1neg;
 128    delete hITSz1;
 129
 130 //-----------------------Get points---------------------------------------------
 131
 132    Double_t dphi,DeltaPhi,DeltaZ,r;
 133    Int_t np1=0;
 134    Int_t np2=0;
 135
 136    Double_t *Z1, *Z2, *Y1, *Y2, *X1, *X2, *phi1, *phi2, *r1, *r2;
 137    Z1=new Double_t[NoPoints1];
 138    Z2=new Double_t[NoPoints2];
 139    Y1=new Double_t[NoPoints1];
 140    Y2=new Double_t[NoPoints2];
 141    X1=new Double_t[NoPoints1];
 142    X2=new Double_t[NoPoints2];
 143    phi1=new Double_t[NoPoints1];
 144    phi2=new Double_t[NoPoints2];
 145    r1=new Double_t[NoPoints1];
 146    r2=new Double_t[NoPoints2];
 147
 148    for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
 149    {
 150         aliits->ResetRecPoints();
 151         gAlice->TreeR()->GetEvent(i);
 152            npoints = recpoints->GetEntries();
 153            for (ipoint=0;ipoint<npoints;ipoint++) {
 154
 155                 //if(rnd.Integer(10)>2) continue;
 156                     pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
 157                 l[0]=pnt->GetX();
 158                 l[1]=0;
 159                 l[2]=pnt->GetZ();
 160                 g2->LtoG(i, l, p);
 161                     p[0]=p[0]-Vzero[0];
 162                     p[1]=p[1]-Vzero[1];
 163                 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 164
 165                         if(i<80 && TMath::Abs(p[2])<14.35) {
 166                         Z1[np1]=p[2];
 167                         Y1[np1]=p[1];
 168                         X1[np1]=p[0];
 169                         r1[np1]=r;
 170                         phi1[np1]=PhiFunc(p);
 171                         np1++;
 172                         }
 173
 174                         if(i>=80 &&  TMath::Abs(p[2])<14.35) {
 175                         Z2[np2]=p[2];
 176                         Y2[np2]=p[1];
 177                         X2[np2]=p[0];
 178                         r2[np2]=r;
 179                         phi2[np2]=PhiFunc(p);
 180                         np2++;
 181                         }
 182
 183           }
 184    }
 185
 186 //-------------------Correlations-----------------------------------------------
 187
 188    //cout << "\nPoints number on the first and second layer: "<<np1<<" "<<np2<<endl;
 189
 190    DeltaPhi = 0.085-0.1*TMath::Sqrt(Vzero[0]*Vzero[0]+Vzero[1]*Vzero[1]);
 191    if(DeltaPhi<0) {
 192       cout << "\nThe algorith can't find the vertex. " << endl;
 193           exit(123456789);
 194    }
 195
 196    Float_t B[3];
 197    Float_t origin[3];
 198    for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
 199    gAlice->Field()->Field(origin,B);
 200
 201    DeltaPhi = DeltaPhi*B[2]/2;
 202
 203    cout << "\nDeltaPhi: " << DeltaPhi << " deg" << "\n";
 204
 205    DeltaZ =
 206    0.2+2.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2);
 207
 208    cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";
 209
 210    TH1F *hITSZv         = new TH1F("hITSZv","",DeltaZ/0.005,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);
 211
 212    for(j=0; j<(np2)-1; j++) {
 213          for(k=0; k<(np1)-1; k++) {
 214                 if(TMath::Abs(Z1[k]-Z2[j])>13.) continue;
 215                 dphi=TMath::Abs(phi1[k]-phi2[j]);
 216                 if(dphi>180) dphi = 360-dphi;
 217                 if(dphi<DeltaPhi) {
 218                   if(TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Vzero[2])<DeltaZ)
 219                       hITSZv->Fill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1));
 220                 }
 221          }
 222    }
 223
 224
 225
 226    cout << "\nNumber of used pairs: \n";
 227    cout << hITSZv->GetEntries() << '\n' << '\n';
 228    cout << "\nCentroid and RMS of Zv distribution: \n";
 229    cout <<  hITSZv->GetMean() << "       " << hITSZv->GetRMS() << "\n"<< "\n";
 230
 231    delete [] Z1;
 232    delete [] Z2;
 233    delete [] Y1;
 234    delete [] Y2;
 235    delete [] X1;
 236    delete [] X2;
 237    delete [] r1;
 238    delete [] r2;
 239    delete [] phi1;
 240    delete [] phi2;
 241
 242    Double_t a = Vzero[2]-DeltaZ;
 243    Double_t b = Vzero[2]+DeltaZ;
 244
 245    TF1 *f4 = new TF1 ("f4","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);
 246    f4->SetParameter(0,700.);
 247    f4->SetParameter(1,Vzero[2]);
 248    f4->SetParameter(2,0.05); // !!
 249    f4->SetParameter(3,50.);
 250
 251    hITSZv->Fit("f4","RME0");
 252
 253
 254    delete hITSZv;
 255
 256    fSNR[2] = f4->GetParameter(0)/f4->GetParameter(3);
 257    if(fSNR[2]<1.5) {
 258      cout << "\nSignal to noise ratio too small!!!" << endl;
 259      cout << "The algorithm can't find the z vertex position." << endl;
 260          exit(123456789);
 261    }
 262    else
 263    {
 264      fPosition[2] = (f4->GetParameter(1));
 265      if(fPosition[2]<0)
 266        {
 267          fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000;
 268        }
 269      else
 270        {
 271          fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000;
 272        }
 273    }
 274    fResolution[2] = f4->GetParError(1);
 275
 276    AliGenerator *gener = gAlice->Generator();
 277    Float_t Vx,Vy,Vz;
 278    gener->GetOrigin(Vx,Vy,Vz);
 279
 280    fPosition[0]=(Double_t)Vx;
 281    fPosition[1]=(Double_t)Vy;
 282
 283    fResolution[0] = 0;
 284    fResolution[1] = 0;
 285
 286    fSNR[0] = 0;
 287    fSNR[1] = 0;
 288
 289 }
 290
 291
 292 //______________________________________________________________________________
 293
 294 AliITSVertex::~AliITSVertex() {
 295    delete [] fPosition;
 296    delete [] fResolution;
 297    delete [] fSNR;
 298 }
 299
 300 //______________________________________________________________________________
 301
 302
 303 Double_t AliITSVertex::PhiFunc(Float_t p[]) {
 304          Double_t phi=0;
 305          if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
 306          if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
 307          if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
 308          if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
 309          return phi;
 310 }
 311 //
 312 //______________________________________________________________________________
 313