ITS/AliITSVertex.cxx

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