ITS/AliITSVertexer3DTapan.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------
  17 //           AliITSVertexer3DTapan class
  18 //   This is a class for the 3d vertex finding
  19 //    Origin: Tapan Nayak, VECC-CERN, Tapan.Nayak@cern.ch
  20 //-----------------------------------------------------------------
  21
  22 #include <TH1.h>
  23 #include <TTree.h>
  24 #include <TClonesArray.h>
  25
  26 #include <AliITSVertexer3DTapan.h>
  27 #include <AliITSRecPoint.h>
  28 #include <AliITSgeom.h>
  29 #include <AliESDVertex.h>
  30
  31 ClassImp(AliITSVertexer3DTapan)
  32
  33 Int_t AliITSVertexer3DTapan::LoadClusters(TTree *cTree) {
  34   //--------------------------------------------------------------------
  35   //This function loads the SPD clusters
  36   //--------------------------------------------------------------------
  37    TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
  38    TBranch *branch=cTree->GetBranch("ITSRecPoints");
  39    branch->SetAddress(&clusters);
  40
  41    Int_t nentr=cTree->GetEntries(),nc1=0,nc2=0;
  42    for (Int_t i=0; i<nentr; i++) {
  43        if (!cTree->GetEvent(i)) continue;
  44        //
  45        //   Below:
  46        //   "alpha" is the angle from the global X-axis to the
  47        //    local GEANT X'-axis  ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
  48        //    "phi" is the angle from the global X-axis to the
  49        //         local cluster X"-axis
  50        //
  51
  52        //       Double_t rot[9];   fITSgeom->GetRotMatrix(i,rot);
  53        Int_t    lay,lad,det; fITSgeom->GetModuleId(i,lay,lad,det);
  54
  55        if (lay>2) break;  //load the SPD clusters only
  56
  57        /*
  58        Float_t  tx,ty,tz;  fITSgeom->GetTrans(lay,lad,det,tx,ty,tz);
  59
  60        Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
  61        Double_t phi=TMath::Pi()/2+alpha;
  62
  63        if (lay==1) phi+=TMath::Pi();
  64        Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
  65        Double_t r=tx*cp+ty*sp;
  66        */
  67
  68        Int_t ncl=clusters->GetEntriesFast();
  69        Float_t hPhi;
  70        while (ncl--) {
  71           AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
  72           Float_t pos[3];
  73           c->GetGlobalXYZ(pos);
  74           if (lay==1) {
  75             /*             fX1[nc1]= r*cp - c->GetY()*sp;
  76              fY1[nc1]= r*sp + c->GetY()*cp;
  77              fZ1[nc1]= c->GetZ(); */
  78             fX1[nc1] = pos[0]; fY1[nc1] = pos[1]; fZ1[nc1] = pos[2];
  79              CalculatePhi(fX1[nc1], fY1[nc1], hPhi);
  80              fPhi1[nc1]= hPhi;
  81              nc1++;
  82           } else {
  83             /* fX2[nc2]= r*cp - c->GetY()*sp;
  84              fY2[nc2]= r*sp + c->GetY()*cp;
  85              fZ2[nc2]= c->GetZ(); */
  86             fX2[nc2] = pos[0]; fY2[nc2] = pos[1]; fZ2[nc2] = pos[2];
  87              CalculatePhi(fX2[nc2], fY2[nc2], hPhi);
  88              fPhi2[nc2]= hPhi;
  89              nc2++;
  90           }
  91        }
  92    }
  93    ficlu1 = nc1; ficlu2 = nc2;
  94    AliInfo(Form("Number of clusters: %d (first layer) and %d (second layer)",ficlu1,ficlu2));
  95    return 0;
  96 }
  97
  98 void AliITSVertexer3DTapan::FindVertexForCurrentEvent(AliESDVertex *vtx) {
  99   //
 100   // This function reconstructs ....
 101   //
 102   //
 103   if (vtx==0) return;
 104
 105   Double_t pos[3], postemp[3], sigpos[3];
 106   Int_t ncontr, ncontrtemp;
 107   Float_t cuts[3];
 108   Int_t vtxstatus=0;
 109
 110   //....
 111   pos[0] = 0.;   pos[1] = 0.;   pos[2] = 0.;
 112   cuts[0]=1.;  cuts[1]=1.;  cuts[2]=20.;
 113   CalculateVertex3d1(pos, cuts, ncontr);
 114   if(ncontr==0) {
 115     pos[0] = 9999.;   pos[1] = 9999.;   pos[2] = 9999.;
 116     vtxstatus = -1;
 117   }
 118   AliInfo(Form("1st step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
 119
 120   if(vtxstatus == 0) {
 121     ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2];
 122     cuts[0]=0.3;  cuts[1]=0.3;  cuts[2]=1.;
 123     CalculateVertex3d1(pos, cuts, ncontr);
 124     if(ncontr==0) {
 125       ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
 126       vtxstatus = 2;
 127     }
 128     AliInfo(Form("2nd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
 129   }
 130
 131   if(vtxstatus == 0) {
 132     ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2];
 133     cuts[0]=0.25;  cuts[1]=0.25;  cuts[2]=1.0;
 134     CalculateVertex3d2(pos, cuts, ncontr, sigpos);
 135     if(ncontr==0) {
 136       ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
 137       vtxstatus = 3;
 138     }
 139     AliInfo(Form("3rd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
 140   }
 141
 142   if(vtxstatus == 0) {
 143     ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2];
 144     cuts[0]=0.1;  cuts[1]=0.1;  cuts[2]=0.2;
 145     CalculateVertex3d2(pos, cuts, ncontr, sigpos);
 146     if(ncontr==0) {
 147       ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
 148       vtxstatus = 4;
 149     }
 150     AliInfo(Form("4th step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
 151   }
 152   AliInfo(Form("Final step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
 153
 154   new(vtx) AliESDVertex(pos,sigpos,(Double_t)vtxstatus,ncontr,"AliITSVertexer3DTapan");
 155   return;
 156 }
 157
 158
 159 void AliITSVertexer3DTapan::CalculateVertex3d1(Double_t pos[3], Float_t cuts[3], Int_t &ncontr) {
 160   //
 161   // This function reconstructs first two steps of vertex
 162   //
 163
 164   Double_t p1[4], p2[4], p3[4], p4[4];
 165   Double_t pa[3], pb[3];
 166   Double_t hphi1, hphi2, hphi3, hphi4;
 167
 168   ncontr = 0;
 169   Float_t  phicut = 1.0;
 170   Double_t distance;  Float_t  distancecut = 1.0;
 171   Int_t    ibin=20;  Float_t  ilow=-1.; Float_t  ihigh=1.;
 172   Int_t    ibinz=400; Float_t  ilowz=-20.; Float_t  ihighz=20.;
 173
 174   TH1F *hx = new TH1F("hx","",  ibin, ilow, ihigh);
 175   TH1F *hy = new TH1F("hy","",  ibin, ilow, ihigh);
 176   TH1F *hz = new TH1F("hz","",  ibinz,ilowz,ihighz);
 177
 178   for (Int_t ip1=0; ip1<ficlu1; ip1++) {
 179     // Two points on layer1: p1 and p3
 180     p1[0] = fX1[ip1];   p1[1] = fY1[ip1];    p1[2] = fZ1[ip1];
 181     p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1];  p3[2] = fZ1[ip1+1];
 182     hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
 183
 184     for (Int_t ip2=0; ip2<ficlu2; ip2++) {
 185       // Two points on layer 2: p2 and p4
 186       p2[0] = fX2[ip2];   p2[1] = fY2[ip2];     p2[2] = fZ2[ip2];
 187       p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1];   p4[2] = fZ2[ip2+1];
 188       hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];
 189
 190       // First line is formed by p1-p2 and second line by p3-p4
 191       // We find two points on each line which form the closest distance of the two lines
 192       // pa[0],pa[1],pa[2]: points on line 1 and pb[0],pb[1],pb[2]: points on line 2
 193       // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
 194
 195       if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
 196         CalculateLine(p1, p2, p3, p4, pa, pb);
 197
 198         if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
 199           distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
 200           if(distance<distancecut){
 201             hx->Fill(pa[0]);            hy->Fill(pa[1]);                hz->Fill(pa[2]);
 202             hx->Fill(pb[0]);            hy->Fill(pb[1]);                hz->Fill(pb[2]);
 203             ncontr++;
 204           }
 205         }
 206       }
 207
 208       // Third line is formed by p1-p4 and fourth line by p3-p2
 209       // We find two points on each line which form the closest distance of the two lines
 210       // pa[0],pa[1],pa[2]: points on line 3 and pb[0],pb[1],pb[2]: points on line 4
 211       // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
 212       if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
 213         CalculateLine(p1, p4, p3, p2, pa, pb);
 214         if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1]){
 215           distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
 216           if(distance<distancecut){
 217             hx->Fill(pa[0]);            hy->Fill(pa[1]);                hz->Fill(pa[2]);
 218             hx->Fill(pb[0]);            hy->Fill(pb[1]);                hz->Fill(pb[2]);
 219             ncontr++;
 220           }
 221         }
 222       }
 223     }
 224   }
 225
 226   Int_t maxbinx = hx->GetMaximumBin();
 227   Int_t maxbiny = hy->GetMaximumBin();
 228   Int_t maxbinz = hz->GetMaximumBin();
 229   pos[0] = ilow  + ((ihigh-ilow)/ibin)*maxbinx;
 230   pos[1] = ilow  + ((ihigh-ilow)/ibin)*maxbiny;
 231   pos[2] = ilowz + ((ihighz-ilowz)/ibinz)*maxbinz;
 232   hx->Delete();
 233   hy->Delete();
 234   hz->Delete();
 235 }
 236
 237 void AliITSVertexer3DTapan::CalculateVertex3d2(Double_t pos[3], Float_t cuts[3], Int_t &ncontr, Double_t sigpos[3]) {
 238   //
 239   // This function reconstructs second two steps of vertex
 240   //
 241
 242   Double_t p1[4], p2[4], p3[4], p4[4];
 243   Double_t pa[3], pb[3];
 244   Double_t hphi1, hphi2, hphi3, hphi4;
 245
 246   ncontr = 0;
 247   Float_t  phicut = 0.3;
 248   Double_t distance; Float_t  distancecut = 1.0;
 249
 250   Double_t vertx  =0.;   Double_t verty  =0.;   Double_t vertz  =0.;
 251   Double_t vertx2 =0.;   Double_t verty2 =0.;   Double_t vertz2 =0.;
 252
 253   for (Int_t ip1=0; ip1<ficlu1; ip1++) {
 254     // Two points on layer1: p1 and p3
 255     p1[0] = fX1[ip1];   p1[1] = fY1[ip1];    p1[2] = fZ1[ip1];
 256     p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1];  p3[2] = fZ1[ip1+1];
 257     hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
 258
 259     for (Int_t ip2=0; ip2<ficlu2; ip2++) {
 260       // Two points on layer 2: p2 and p4
 261       p2[0] = fX2[ip2];   p2[1] = fY2[ip2];     p2[2] = fZ2[ip2];
 262       p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1];   p4[2] = fZ2[ip2+1];
 263       hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];
 264
 265       // First line is formed by p1-p2 and second line by p3-p4
 266       // We find two points on each line which form the closest distance of the two lines
 267       // pa[0],pa[1],pa[2] are the points on line 1 and pb[0],pb[1],pb[2] are the points on line 2
 268       // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
 269
 270       if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
 271         CalculateLine(p1, p2, p3, p4, pa, pb);
 272
 273         // We consider the points where x, y and z points are less than xcut, ycut and zcut, respectively
 274         if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
 275           distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
 276           if(distance<distancecut){
 277             ncontr++;
 278             vertx  = vertx  + pa[0];       verty  = verty  + pa[1];       vertz  = vertz  + pa[2];
 279             vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
 280             ncontr++;
 281             vertx  = vertx  + pb[0];       verty  = verty  + pb[1];       vertz  = vertz  + pb[2];
 282             vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
 283           }
 284         }
 285       }
 286
 287       // Third line is formed by p1-p4 and fourth line by p3-p2
 288       // We find two points on each line which form the closest distance of the two lines
 289       // pa[0],pa[1],pa[2] are the points on line 3 and pb[0],pb[1],pb[2] are the points on line 4
 290       // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
 291       if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
 292
 293         CalculateLine(p1, p4, p3, p2, pa, pb);
 294         if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
 295           distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
 296           if(distance<distancecut){
 297             ncontr++;
 298             vertx  = vertx  + pa[0];       verty  = verty  + pa[1];       vertz  = vertz  + pa[2];
 299             vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
 300             ncontr++;
 301             vertx  = vertx  + pb[0];       verty  = verty  + pb[1];       vertz  = vertz  + pb[2];
 302             vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
 303           }
 304         }
 305       }
 306     }
 307   }
 308
 309   if(ncontr>0){
 310     pos[0]  = vertx/ncontr;  pos[1] = verty/ncontr;   pos[2] = vertz/ncontr;
 311     vertx2 = vertx2/ncontr; verty2 = verty2/ncontr; vertz2 = vertz2/ncontr;
 312     sigpos[0] = TMath::Sqrt(vertx2 - pos[0]*pos[0]);
 313     sigpos[1] = TMath::Sqrt(verty2 - pos[1]*pos[1]);
 314     sigpos[2] = TMath::Sqrt(vertz2 - pos[2]*pos[2]);
 315   }
 316   ncontr = ncontr/2;
 317 }
 318
 319 void AliITSVertexer3DTapan::CalculatePhi(Float_t fx, Float_t fy, Float_t & phi)
 320 {
 321   //calculates phi
 322   Float_t ybyx, phi1;
 323   const Float_t kradian = 180./3.141592654;
 324
 325   if(fx==0.)
 326     {
 327       if(fy>0.) phi = 90.;
 328       if(fy<0.) phi = 270.;
 329     }
 330   if(fx != 0.)
 331     {
 332       ybyx = fy/fx;
 333       if(ybyx < 0) ybyx = - ybyx;
 334       phi1 = TMath::ATan(ybyx)*kradian;
 335       if(fx > 0 && fy > 0) phi = phi1;        // 1st Quadrant
 336       if(fx < 0 && fy > 0) phi = 180 - phi1;  // 2nd Quadrant
 337       if(fx < 0 && fy < 0) phi = 180 + phi1;  // 3rd Quadrant
 338       if(fx > 0 && fy < 0) phi = 360 - phi1;  // 4th Quadrant
 339
 340     }
 341   phi = phi/kradian;
 342 }
 343
 344 void AliITSVertexer3DTapan::CalculateLine(Double_t p1[4], Double_t p2[4], Double_t p3[4], Double_t p4[4], Double_t pa[3], Double_t pb[3]) const{
 345   //calculates line
 346   Double_t p13x, p13y, p13z;
 347   Double_t p21x, p21y, p21z;
 348   Double_t p43x, p43y, p43z;
 349   Double_t d1343, d4321, d1321, d4343, d2121;
 350   Double_t numer, denom;
 351   Double_t mua,   mub;
 352   mua = 0.; mub = 0.;
 353
 354   p13x = p1[0] - p3[0];
 355   p13y = p1[1] - p3[1];
 356   p13z = p1[2] - p3[2];
 357
 358   p21x = p2[0] - p1[0];
 359   p21y = p2[1] - p1[1];
 360   p21z = p2[2] - p1[2];
 361
 362   p43x = p4[0] - p3[0];
 363   p43y = p4[1] - p3[1];
 364   p43z = p4[2] - p3[2];
 365
 366   d1343 = p13x * p43x + p13y * p43y + p13z * p43z;
 367   d4321 = p43x * p21x + p43y * p21y + p43z * p21z;
 368   d1321 = p13x * p21x + p13y * p21y + p13z * p21z;
 369   d4343 = p43x * p43x + p43y * p43y + p43z * p43z;
 370   d2121 = p21x * p21x + p21y * p21y + p21z * p21z;
 371
 372   denom = d2121 * d4343 - d4321 * d4321;
 373   numer = d1343 * d4321 - d1321 * d4343;
 374
 375   if(denom>0) mua = numer / denom;
 376   if(d4343>0) mub = (d1343 + d4321 * (mua)) / d4343;
 377
 378   pa[0] = p1[0] + mua * p21x;
 379   pa[1] = p1[1] + mua * p21y;
 380   pa[2] = p1[2] + mua * p21z;
 381
 382   pb[0] = p3[0] + mub * p43x;
 383   pb[1] = p3[1] + mub * p43y;
 384   pb[2] = p3[2] + mub * p43z;
 385 }
 386