RICH/AliRICHPatRec.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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   $Log$
  18   Revision 1.12  2001/05/10 12:34:23  jbarbosa
  19   Changed drwaing routines.
  20
  21   Revision 1.11  2001/03/14 18:21:24  jbarbosa
  22   Corrected bug in digits loading.
  23
  24   Revision 1.10  2001/02/27 15:21:06  jbarbosa
  25   Transition to SDigits.
  26
  27   Revision 1.9  2001/02/13 20:38:48  jbarbosa
  28   Changes to make it work with new IO.
  29
  30   Revision 1.8  2000/11/01 15:37:18  jbarbosa
  31   Updated to use its own rec. point object.
  32
  33   Revision 1.7  2000/10/03 21:44:09  morsch
  34   Use AliSegmentation and AliHit abstract base classes.
  35
  36   Revision 1.6  2000/10/02 21:28:12  fca
  37   Removal of useless dependecies via forward declarations
  38
  39   Revision 1.5  2000/10/02 15:50:25  jbarbosa
  40   Fixed forward declarations.
  41
  42   Revision 1.4  2000/06/30 16:33:43  dibari
  43   Several changes (ring drawing, fiducial selection, etc.)
  44
  45   Revision 1.3  2000/06/15 15:47:12  jbarbosa
  46   Corrected compilation errors on HP-UX (replaced pow with TMath::Power)
  47
  48   Revision 1.2  2000/06/12 15:26:09  jbarbosa
  49   Cleaned up version.
  50
  51   Revision 1.1  2000/06/09 14:53:01  jbarbosa
  52   Bari's pattern recognition algorithm
  53
  54 */
  55
  56 #include "AliRICHHit.h"
  57 #include "AliRICHCerenkov.h"
  58 #include "AliRICHSDigit.h"
  59 #include "AliRICHDigit.h"
  60 #include "AliRICHRawCluster.h"
  61 #include "AliRICHRecHit1D.h"
  62 #include "AliRun.h"
  63 #include "AliDetector.h"
  64 #include "AliRICH.h"
  65 #include "AliRICHPoints.h"
  66 #include "AliSegmentation.h"
  67 #include "AliRICHPatRec.h"
  68 #include "AliRICH.h"
  69 #include "AliRICHConst.h"
  70 #include "AliRICHPoints.h"
  71 #include "AliConst.h"
  72 #include "AliHitMap.h"
  73
  74 #include <TParticle.h>
  75 #include <TMath.h>
  76 #include <TRandom.h>
  77 #include <TCanvas.h>
  78 #include <TH2.h>
  79 #include <TTree.h>
  80
  81
  82 ClassImp(AliRICHPatRec)
  83 //___________________________________________
  84 AliRICHPatRec::AliRICHPatRec() : TObject()
  85 {
  86   // Default constructor
  87
  88     //fChambers = 0;
  89 }
  90 //___________________________________________
  91 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
  92     : TObject()
  93 {
  94     //Constructor for Bari's pattern recogniton method object
  95 }
  96
  97 void AliRICHPatRec::PatRec()
  98 {
  99
 100 // Pattern recognition algorithm
 101
 102   AliRICHChamber*       iChamber;
 103   AliSegmentation*  segmentation;
 104
 105   Int_t ntracks, ndigits[kNCH];
 106   Int_t itr, ich, i;
 107   Int_t goodPhotons;
 108   Int_t x,y,q;
 109   Float_t rx,ry,rz;
 110   Int_t nent,status;
 111   Int_t padsUsedX[100];
 112   Int_t padsUsedY[100];
 113
 114   Float_t rechit[7];
 115
 116   //printf("PatRec started\n");
 117
 118   AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
 119   TTree *treeH = gAlice->TreeH();
 120
 121   ntracks =(Int_t) treeH->GetEntries();
 122   //  ntracks = 1;
 123   for (itr=0; itr<ntracks; itr++) {
 124
 125     status = TrackParam(itr,ich,0,0);
 126     if(status==1) continue;
 127     //printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
 128     //    ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
 129
 130     iChamber = &(pRICH->Chamber(ich));
 131     segmentation=iChamber->GetSegmentationModel();
 132
 133     nent=(Int_t)gAlice->TreeD()->GetEntries();
 134     gAlice->TreeD()->GetEvent(0);
 135     TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
 136     ndigits[ich] = pDigitss->GetEntriesFast();
 137     printf("Digits in chamber %d: %d\n",ich,ndigits[ich]);
 138     AliRICHDigit *padI = 0;
 139
 140     goodPhotons = 0;
 141
 142     for (Int_t dig=0;dig<ndigits[ich];dig++) {
 143       padI=(AliRICHDigit*) pDigitss->UncheckedAt(dig);
 144       x=padI->PadX();
 145       y=padI->PadY();
 146       q=padI->Signal();
 147       segmentation->GetPadC(x,y,rx,ry,rz);
 148
 149       //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
 150       //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
 151
 152       fXpad = rx-fXshift;
 153       fYpad = ry-fYshift;
 154       fQpad = q;
 155
 156       fCerenkovAnglePad = PhotonCerenkovAngle();
 157       if(fCerenkovAnglePad==-999) continue;
 158
 159       if(!PhotonInBand()) continue;
 160
 161       Int_t xpad;
 162       Int_t ypad;
 163
 164       segmentation->GetPadI(fXpad,fYpad,0,xpad,ypad);
 165
 166       padsUsedX[goodPhotons]=xpad;
 167       padsUsedY[goodPhotons]=ypad;
 168
 169       goodPhotons++;
 170       fEtaPhotons[goodPhotons-1] = fCerenkovAnglePad;
 171     }
 172     fNumEtaPhotons = goodPhotons;
 173
 174     BackgroundEstimation();
 175
 176     HoughResponse();
 177     //CerenkovRingDrawing();
 178
 179     rechit[0] = 0;
 180     rechit[1]   = 0;
 181     rechit[2] = fThetaCerenkov;
 182     rechit[3] = fXshift + fTrackLoc[0];
 183     rechit[4] = fYshift + fTrackLoc[1];
 184     rechit[5] = fEmissPoint;
 185     rechit[6] = goodPhotons;
 186
 187     //printf("Center coordinates:%f %f\n",rechit[3],rechit[4]);
 188
 189     pRICH->AddRecHit1D(ich,rechit,fEtaPhotons,padsUsedX,padsUsedY);
 190
 191   }
 192
 193   gAlice->TreeR()->Fill();
 194   TClonesArray *fRec;
 195   for (i=0;i<kNCH;i++) {
 196     fRec=pRICH->RecHitsAddress1D(i);
 197     int ndig=fRec->GetEntriesFast();
 198     printf ("Chamber %d, rings %d\n",i,ndig);
 199   }
 200   pRICH->ResetRecHits1D();
 201
 202 }
 203
 204
 205 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich, Float_t rectheta, Float_t recphi)
 206 {
 207   // Get Local coordinates of track impact
 208
 209   AliRICHChamber*       iChamber;
 210   AliSegmentation*      segmentation;
 211
 212   Float_t trackglob[3];
 213   Float_t trackloc[3];
 214   Float_t thetatr;
 215   Float_t phitr;
 216   Float_t iloss;
 217   Float_t part;
 218   Float_t pX, pY, pZ;
 219
 220   //printf("Calling TrackParam\n");
 221
 222     gAlice->ResetHits();
 223     TTree *treeH = gAlice->TreeH();
 224     treeH->GetEvent(itr);
 225
 226     AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
 227     AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
 228     if(mHit==0) return 1;
 229     ich = mHit->Chamber()-1;
 230     trackglob[0] = mHit->X();
 231     trackglob[1] = mHit->Y();
 232     trackglob[2] = mHit->Z();
 233     pX = mHit->MomX();
 234     pY = mHit->MomY();
 235     pZ = mHit->MomZ();
 236     fTrackMom = sqrt(TMath::Power(pX,2)+TMath::Power(pY,2)+TMath::Power(pZ,2));
 237     if(recphi!=0 || rectheta!=0)
 238       {
 239         thetatr = rectheta;
 240         phitr = recphi;
 241       }
 242     else
 243       {
 244         thetatr = mHit->Theta()*TMath::Pi()/180;
 245         phitr = mHit->Phi()*TMath::Pi()/180;
 246       }
 247     iloss = mHit->Loss();
 248     part  = mHit->Particle();
 249
 250     iChamber = &(pRICH->Chamber(ich));
 251     iChamber->GlobaltoLocal(trackglob,trackloc);
 252
 253     segmentation=iChamber->GetSegmentationModel();
 254
 255     // retrieve geometrical params
 256
 257     AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();
 258
 259     fRw   = fGeometry->GetFreonThickness();
 260     fQw   = fGeometry->GetQuartzThickness();
 261     fTgap = fGeometry->GetGapThickness();
 262     Float_t radiatorToPads= fGeometry->GetRadiatorToPads();
 263       //+ fGeometry->GetProximityGapThickness();
 264
 265     //printf("Distance to pads. From geometry:%f, From calculations:%f\n",radiatorToPads,fRw + fQw + fTgap);
 266
 267     //Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
 268     Float_t apar = radiatorToPads*tan(thetatr);
 269     fTrackLoc[0] = apar*cos(phitr);
 270     fTrackLoc[1] = apar*sin(phitr);
 271     //fTrackLoc[2] = fRw + fQw + fTgap;
 272     fTrackLoc[2] = radiatorToPads;
 273     fTrackTheta  = thetatr;
 274     fTrackPhi    = phitr;
 275
 276     fXshift = trackloc[0] - fTrackLoc[0];
 277     fYshift = trackloc[2] - fTrackLoc[1];
 278
 279     return 0;
 280 }
 281
 282 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
 283                                                  Float_t phiphot)
 284 {
 285
 286 // Estimation of emission point
 287
 288   Float_t nquartz = 1.585;
 289   Float_t ngas    = 1.;
 290   Float_t nfreon  = 1.295;
 291   Float_t value;
 292
 293   //  printf("Calling EstimationLimits\n");
 294
 295   Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
 296   Float_t b1 = (fRw-fEmissPoint)*tan(lim);
 297   Float_t b2 = fQw / sqrt(TMath::Power(nquartz,2)-TMath::Power(nfreon*sin(lim),2));
 298   Float_t b3 = fTgap / sqrt(TMath::Power(ngas,2)-TMath::Power(nfreon*sin(lim),2));
 299   Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
 300   value = TMath::Power(radius,2)
 301     -TMath::Power((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
 302     -TMath::Power((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
 303   return value;
 304 }
 305
 306
 307 Float_t AliRICHPatRec::PhotonCerenkovAngle()
 308 {
 309   // Cherenkov pad angle reconstruction
 310
 311   Float_t radius;
 312   Float_t cherMin = 0;
 313   Float_t cherMax = 0.8;
 314   Float_t phiphot;
 315   Float_t eps = 0.0001;
 316   Int_t   niterEmiss = 0;
 317   Int_t   niterEmissMax = 0;
 318   Float_t x1,x2,x3=0,p1,p2,p3;
 319   Float_t argY,argX;
 320   Int_t niterFun;
 321
 322   //  printf("Calling PhotonCerenkovAngle\n");
 323
 324   radius = sqrt(TMath::Power(fTrackLoc[0]-fXpad,2)+TMath::Power(fTrackLoc[1]-fYpad,2));
 325   fEmissPoint = fRw/2.;  //Start value of EmissionPoint
 326
 327   while(niterEmiss<=niterEmissMax) {
 328
 329     niterFun = 0;
 330     argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
 331     argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
 332     phiphot = atan2(argY,argX);
 333     p1 = EstimationAtLimits(cherMin,radius,phiphot);
 334     p2 = EstimationAtLimits(cherMax,radius,phiphot);
 335     if(p1*p2>0)
 336       {
 337         //        printf("PhotonCerenkovAngle failed\n");
 338         return -999;
 339       }
 340
 341     //start to find the Cherenkov pad angle
 342     x1 = cherMin;
 343     x2 = cherMax;
 344     x3 = (x1+x2)/2.;
 345     p3 = EstimationAtLimits(x3,radius,phiphot);
 346     while(TMath::Abs(p3)>eps){
 347       if(p1*p3<0) x2 = x3;
 348       if(p1*p3>0) {
 349         x1 = x3;
 350         p1 = EstimationAtLimits(x1,radius,phiphot);
 351       }
 352       x3 = (x1+x2)/2.;
 353       p3 = EstimationAtLimits(x3,radius,phiphot);
 354       niterFun++;
 355
 356       if(niterFun>=1000) {
 357         //      printf(" max iterations in PhotonCerenkovAngle\n");
 358         return x3;
 359       }
 360     }
 361     //    printf("niterFun %i \n",niterFun);
 362     niterEmiss++;
 363     if (niterEmiss != niterEmissMax+1) EmissionPoint();
 364   }
 365   /*
 366    printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
 367                          phiphot,fXpad,fYpad,fEmissPoint);
 368   */
 369
 370   return x3;
 371
 372 }
 373
 374
 375 void AliRICHPatRec::EmissionPoint()
 376 {
 377
 378 // Find emission point
 379
 380   Float_t absorbtionLength=7.83*fRw; //absorption length in the freon (cm)
 381   // 7.83 = -1/ln(T0) where
 382   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
 383   Float_t photonLength, photonLengthMin, photonLengthMax;
 384
 385   photonLength=exp(-fRw/(absorbtionLength*cos(fCerenkovAnglePad)));
 386   photonLengthMin=fRw*photonLength/(1.-photonLength);
 387   photonLengthMax=absorbtionLength*cos(fCerenkovAnglePad);
 388   fEmissPoint = fRw + photonLengthMin - photonLengthMax;
 389
 390 }
 391
 392 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
 393 {
 394
 395 // not implemented yet
 396
 397   printf("Calling PhotonSelection\n");
 398 }
 399
 400 void AliRICHPatRec::BackgroundEstimation()
 401 {
 402
 403 // estimate background noise
 404
 405   Float_t stepEta   = 0.001;
 406   Float_t etaMinBkg = 0.72;
 407   Float_t etaMaxBkg = 0.75;
 408   Float_t etaMin    = 0.;
 409   Float_t etaMax    = 0.75;
 410   Float_t ngas      = 1.;
 411   Float_t nfreon    = 1.295;
 412
 413   Float_t etaStepMin,etaStepMax,etaStepAvg;
 414   Int_t i,ip,nstep;
 415   Int_t numPhotBkg, numPhotonStep;
 416   Float_t funBkg,areaBkg,normBkg;
 417   Float_t densityBkg,storeBkg,numStore;
 418   Float_t thetaSig;
 419
 420   numPhotBkg = 0;
 421   areaBkg = 0.;
 422
 423   nstep = (int)((etaMaxBkg-etaMinBkg)/stepEta);
 424
 425   for (i=0;i<fNumEtaPhotons;i++) {
 426
 427     if(fEtaPhotons[i]>etaMinBkg && fEtaPhotons[i]<etaMaxBkg) {
 428       numPhotBkg++;
 429     }
 430   }
 431   if (numPhotBkg == 0) {
 432      for (i=0;i<fNumEtaPhotons;i++) {
 433         fWeightPhotons[i] = 1.;
 434      }
 435     return;
 436   }
 437
 438   //  printf(" numPhotBkg %i ",numPhotBkg);
 439
 440   for (i=0;i<nstep;i++) {
 441     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
 442     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
 443     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
 444     /*
 445     funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
 446                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
 447     */
 448
 449     //printf("etaStepAvg: %f, etaStepMax: %f, etaStepMin: %f", etaStepAvg,etaStepMax,etaStepMin);
 450
 451     thetaSig = TMath::ASin(nfreon/ngas*TMath::Sin(etaStepAvg));
 452     funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
 453        /ngas*cos(etaStepAvg)/cos(thetaSig);
 454     areaBkg += stepEta*funBkg;
 455   }
 456
 457   densityBkg = 0.95*(Float_t)(numPhotBkg)/areaBkg;
 458   //  printf(" densityBkg %f \n",densityBkg);
 459
 460   nstep = (int)((etaMax-etaMin)/stepEta);
 461   storeBkg = 0.;
 462   numStore = 0;
 463   for (i=0;i<nstep;i++) {
 464     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
 465     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
 466     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
 467     /*
 468     funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
 469                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
 470     */
 471
 472     thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
 473     funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
 474        /ngas*cos(etaStepAvg)/cos(thetaSig);
 475
 476     areaBkg = stepEta*funBkg;
 477     normBkg = densityBkg*areaBkg;
 478     numPhotonStep = 0;
 479     for (ip=0;ip<fNumEtaPhotons;ip++) {
 480       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
 481         numPhotonStep++;
 482       }
 483     }
 484     if (numPhotonStep == 0) {
 485       storeBkg += normBkg;
 486       numStore++;
 487       if (numStore>50) {
 488         numStore = 0;
 489         storeBkg = 0.;
 490       }
 491     }
 492     if (numPhotonStep == 0) continue;
 493     for (ip=0;ip<fNumEtaPhotons;ip++) {
 494       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
 495         normBkg +=storeBkg;
 496         storeBkg = 0;
 497         numStore = 0;
 498         fWeightPhotons[ip] = 1. - normBkg/(Float_t)(numPhotonStep);
 499         /*
 500         printf(" normBkg %f numPhotonStep %i fW %f \n",
 501                normBkg, numPhotonStep, fWeightPhotons[ip]);
 502         */
 503         if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
 504       }
 505     }
 506   }
 507 }
 508
 509
 510 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
 511 {
 512
 513 // not implemented yet
 514
 515   printf("Calling FlagPhotons\n");
 516 }
 517
 518
 519 //////////////////////////////////////////
 520
 521
 522
 523
 524
 525 Int_t AliRICHPatRec::PhotonInBand()
 526 {
 527   //0=label for parameters giving internal band ellipse
 528   //1=label for parameters giving external band ellipse
 529
 530   Float_t imp[2], mass[2], energy[2], beta[2];
 531   Float_t emissPointLength[2];
 532   Float_t e1, e2, f1, f2;
 533   Float_t nfreon[2], nquartz[2];
 534   Int_t times;
 535   Float_t pointsOnCathode[3];
 536
 537   Float_t phpad, thetacer[2];
 538   Float_t bandradius[2], padradius;
 539
 540   imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
 541   imp[1] = 1.2;
 542
 543   mass[0] = 0.938; //proton mass
 544   mass[1] = 0.139; //pion mass
 545
 546   emissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
 547   emissPointLength[1] = 0.;//at the end of radiator
 548
 549   //parameters to calculate freon window refractive index vs. energy
 550   Float_t a = 1.177;
 551   Float_t b = 0.0172;
 552
 553   //parameters to calculate quartz window refractive index vs. energy
 554   /*
 555   Energ[0]  = 5.6;
 556   Energ[1]  = 7.7;
 557   */
 558   energy[0]  = 5.0;
 559   energy[1]  = 8.0;
 560   e1 = 10.666;
 561   e2  = 18.125;
 562   f1  = 46.411;
 563   f2  = 228.71;
 564
 565   phpad = PhiPad(fTrackTheta,fTrackPhi);
 566
 567   for (times=0; times<=1; times++) {
 568
 569     nfreon[times]   = a+b*energy[times];
 570     //nfreon[times]   = 1;
 571
 572     nquartz[times] = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(energy[times],2)))+
 573                           (f2/(TMath::Power(e2,2)-TMath::Power(energy[times],2))));
 574
 575     beta[times]  = imp[times]/sqrt(TMath::Power(imp[times],2)+TMath::Power(mass[times],2));
 576
 577     thetacer[times] =  CherenkovAngle( nfreon[times], beta[times]);
 578
 579     bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
 580                                 emissPointLength[times],
 581                                       thetacer[times], phpad, pointsOnCathode,fTrackTheta,fTrackPhi);
 582     //printf(" ppp %f %f %f \n",pointsOnCathode);
 583 }
 584
 585   bandradius[0] -= 1.6;
 586   bandradius[1] += 1.6;
 587   padradius = sqrt(TMath::Power(fXpad,2)+TMath::Power(fYpad,2));
 588   //  printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
 589
 590   if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
 591   return 0;
 592 }
 593
 594 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
 595                        Float_t emissPointLength, Float_t thetacer,
 596                        Float_t phpad, Float_t pointsOnCathode[3], Float_t rectheta, Float_t recphi)
 597 {
 598
 599 // Find the distance to MIP impact
 600
 601   Float_t distanceValue;
 602
 603   TVector3 radExitPhot(1,1,1);//photon impact at the radiator exit with respect
 604   //to local reference sistem with the origin in the MIP entrance
 605
 606   TVector3 vectEmissPointLength(1,1,1);
 607   Float_t magEmissPointLenght;
 608
 609   TVector3 radExitPhot2(1,1,1);//photon impact at the radiator exit with respect
 610   Float_t magRadExitPhot2;
 611   //to a reference sistem with origin in the photon emission point and
 612   //axes parallel to the MIP reference sistem
 613
 614   TVector3 quarExitPhot(1,1,1);//photon impact at the quartz exit with respect
 615   Float_t magQuarExitPhot;
 616   //
 617   TVector3 gapExitPhot(1,1,1) ;
 618   Float_t magGapExitPhot;
 619   //
 620   TVector3 PhotocatExitPhot(1,1,1);
 621   Double_t theta2;
 622   Double_t thetarad , phirad ;
 623   Double_t thetaquar, phiquar;
 624   Double_t thetagap , phigap ;
 625
 626   Float_t ngas    = 1.;
 627
 628   magEmissPointLenght =  emissPointLength/cos(rectheta);
 629
 630   vectEmissPointLength.SetMag(magEmissPointLenght);
 631   vectEmissPointLength.SetTheta(rectheta);
 632   vectEmissPointLength.SetPhi(recphi);
 633
 634
 635   radExitPhot2.SetTheta(thetacer);
 636   radExitPhot2.SetPhi(phpad);
 637
 638
 639   TRotation r1;
 640   TRotation r2;
 641   TRotation r;
 642
 643   r1. RotateY(rectheta);
 644   r2. RotateZ(recphi);
 645
 646
 647
 648   r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
 649   //following by a rotation about the z axis by MIP phi incidence angle;
 650
 651
 652   radExitPhot2    = r * radExitPhot2;
 653   theta2          = radExitPhot2.Theta();
 654   magRadExitPhot2 = (fRw -  vectEmissPointLength(2))/cos(theta2);
 655   radExitPhot2.SetMag(magRadExitPhot2);
 656
 657
 658   radExitPhot = vectEmissPointLength + radExitPhot2;
 659   thetarad    = radExitPhot.Theta();
 660   phirad  =  radExitPhot.Phi(); //check on the original file //
 661
 662   thetaquar   = SnellAngle( nfreon, nquartz, theta2);
 663   phiquar     = radExitPhot2.Phi();
 664   if(thetaquar == 999.) return thetaquar;
 665   magQuarExitPhot    = fQw/cos(thetaquar);
 666   quarExitPhot.SetMag( magQuarExitPhot);
 667   quarExitPhot.SetTheta(thetaquar);
 668   quarExitPhot.SetPhi(phiquar);
 669
 670   thetagap = SnellAngle( nquartz, ngas, thetaquar);
 671   phigap   = phiquar;
 672   if(thetagap == 999.) return thetagap;
 673   magGapExitPhot    = fTgap/cos(thetagap);
 674   gapExitPhot.SetMag( magGapExitPhot);
 675   gapExitPhot.SetTheta(thetagap);
 676   gapExitPhot.SetPhi(phigap);
 677
 678   PhotocatExitPhot =  radExitPhot + quarExitPhot + gapExitPhot;
 679
 680   distanceValue = sqrt(TMath::Power(PhotocatExitPhot(0),2)
 681                            +TMath::Power(PhotocatExitPhot(1),2));
 682   pointsOnCathode[0] = (Float_t) PhotocatExitPhot(0) + fXshift - fTrackLoc[0];
 683   pointsOnCathode[1] = (Float_t) PhotocatExitPhot(1) + fYshift - fTrackLoc[1];
 684   pointsOnCathode[2] = (Float_t) PhotocatExitPhot(2);
 685
 686   //printf(" point in Distance.2. %f %f %f \n",pointsOnCathode[0],pointsOnCathode[1],pointsOnCathode[2]);
 687
 688   return distanceValue;
 689
 690  }
 691
 692 Float_t AliRICHPatRec::PhiPad(Float_t rectheta, Float_t recphi)
 693 {
 694
 695 // ??
 696
 697   Float_t zpad;
 698   Float_t thetapad, phipad;
 699   Float_t thetarot, phirot;
 700
 701   zpad = fRw + fQw + fTgap;
 702
 703   TVector3 photonPad(fXpad, fYpad, zpad);
 704   thetapad = photonPad.Theta();
 705   phipad = photonPad.Phi();
 706
 707   TRotation r1;
 708   TRotation r2;
 709   TRotation r;
 710
 711   thetarot = - rectheta;
 712   phirot   = - recphi;
 713   r1. RotateZ(phirot);
 714   r2. RotateY(thetarot);
 715
 716   r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
 717   //following by a rotation about the y axis by MIP -theta incidence angle;
 718
 719   photonPad  = r * photonPad;
 720
 721   phipad = photonPad.Phi();
 722
 723   return phipad;
 724 }
 725
 726 Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
 727 {
 728
 729 // Compute the Snell angle
 730
 731   Float_t sinrefractangle;
 732   Float_t refractangle;
 733
 734   sinrefractangle = (n1/n2)*sin(theta1);
 735
 736   if(sinrefractangle>1.) {
 737      refractangle = 999.;
 738      return refractangle;
 739   }
 740
 741   refractangle = asin(sinrefractangle);
 742   return refractangle;
 743 }
 744
 745 Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
 746 {
 747
 748 // Compute the cerenkov angle
 749
 750   Float_t thetacer;
 751
 752   if((n*beta)<1.) {
 753     thetacer = 999.;
 754     return thetacer;
 755   }
 756
 757   thetacer = acos (1./(n*beta));
 758   return thetacer;
 759 }
 760
 761 Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
 762 {
 763
 764 // Find beta
 765
 766   Float_t beta;
 767
 768   beta = 1./(n*cos(theta));
 769   return beta;
 770 }
 771
 772
 773
 774
 775 void AliRICHPatRec::HoughResponse()
 776
 777 {
 778
 779 // Implement Hough response pat. rec. method
 780
 781   int           bin=0;
 782   int           bin1=0;
 783   int           bin2=0;
 784   int           i, j, k, nCorrBand;
 785   int           etaBin = 750;
 786   float         hcs[750];
 787   float         angle, thetaCerMean;
 788
 789   float         etaPeak[30];
 790   float         etaMin = 0.00;
 791   float         etaMax = 0.75;
 792   float         stepEta = 0.001;
 793   float         windowEta = 0.040;
 794
 795   int           nBin;
 796
 797   float etaPeakPos  = -1;
 798   Int_t   etaPeakCount = -1;
 799
 800   thetaCerMean   = 0.;
 801   fThetaCerenkov = 0.;
 802
 803   nBin = (int)(0.5+etaMax/(stepEta));
 804   nCorrBand = (int)(0.5+ windowEta/(2 * stepEta));
 805   memset ((void *)hcs, 0, etaBin*sizeof(int));
 806
 807   for (k=0; k< fNumEtaPhotons; k++) {
 808
 809     angle = fEtaPhotons[k];
 810
 811     if (angle>=etaMin && angle<= etaMax) {
 812       bin = (int)(0.5+angle/(stepEta));
 813       bin1= bin-nCorrBand;
 814       bin2= bin+nCorrBand;
 815       if (bin1<0)    bin1=0;
 816       if (bin2>nBin) bin2=nBin;
 817
 818       for (j=bin1; j<bin2; j++) {
 819         hcs[j] += fWeightPhotons[k];
 820       }
 821
 822       thetaCerMean += angle;
 823     }
 824   }
 825
 826  thetaCerMean /= fNumEtaPhotons;
 827
 828   HoughFiltering(hcs);
 829
 830   for (bin=0; bin <nBin; bin++) {
 831     angle = (bin+0.5) * (stepEta);
 832     if (hcs[bin] && hcs[bin] > etaPeakPos) {
 833       etaPeakCount = 0;
 834       etaPeakPos = hcs[bin];
 835       etaPeak[0]=angle;
 836     }
 837     else {
 838       if (hcs[bin] == etaPeakPos) {
 839         etaPeak[++etaPeakCount] = angle;
 840       }
 841     }
 842   }
 843
 844   for (i=0; i<etaPeakCount+1; i++) {
 845     fThetaCerenkov += etaPeak[i];
 846   }
 847   if (etaPeakCount>=0) {
 848     fThetaCerenkov /= etaPeakCount+1;
 849     fThetaPeakPos = etaPeakPos;
 850   }
 851 }
 852
 853
 854 void AliRICHPatRec::HoughFiltering(float hcs[])
 855 {
 856
 857 // hough filtering
 858
 859    float hcsFilt[750];
 860    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
 861    int nx, i, nxDx;
 862    int sizeHCS;
 863    int nBin;
 864
 865    int   etaBin = 750;
 866    float etaMax = 0.75;
 867    float stepEta = 0.001;
 868
 869    nBin =  (int)(1+etaMax/stepEta);
 870    sizeHCS = etaBin*sizeof(float);
 871
 872    memset ((void *)hcsFilt, 0, sizeHCS);
 873
 874    for (nx = 0; nx < nBin; nx++) {
 875       for (i = 0; i < 5; i++)   {
 876         nxDx = nx + (i-2);
 877         if (nxDx> -1 && nxDx<nBin)
 878              hcsFilt[nx] +=  hcs[nxDx] * k[i];
 879       }
 880    }
 881
 882    for (nx = 0; nx < nBin; nx++) {
 883      hcs[nx] = hcsFilt[nx];
 884    }
 885 }
 886
 887 /*void AliRICHPatRec::CerenkovRingDrawing()
 888 {
 889
 890 //to draw Cherenkov ring by known Cherenkov angle
 891
 892     Int_t nmaxdegrees;
 893     Int_t Nphpad;
 894     Float_t phpad;
 895     Float_t nfreonave, nquartzave;
 896     Float_t aveEnerg;
 897     Float_t energy[2];
 898     Float_t e1, e2, f1, f2;
 899     Float_t bandradius;
 900
 901 //parameters to calculate freon window refractive index vs. energy
 902
 903     Float_t a = 1.177;
 904     Float_t b = 0.0172;
 905
 906 //parameters to calculate quartz window refractive index vs. energy
 907
 908     energy[0]  = 5.0;
 909     energy[1]  = 8.0;
 910     e1  = 10.666;
 911     e2  = 18.125;
 912     f1  = 46.411;
 913     f2  = 228.71;
 914
 915
 916    nmaxdegrees = 36;
 917
 918    for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) {
 919
 920        phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
 921
 922        aveEnerg =  (energy[0]+energy[1])/2.;
 923
 924        nfreonave  = a+b*aveEnerg;
 925        nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
 926                          (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
 927
 928        bandradius = DistanceFromMip(nfreonave, nquartzave,
 929                                    fEmissPoint,fThetaCerenkov, phpad);
 930
 931        fCoordEllipse[0][Nphpad] = fOnCathode[0];
 932        fCoordEllipse[1][Nphpad] = fOnCathode[1];
 933        printf(" values %f %f \n",fOnCathode[0],fOnCathode[1]);
 934
 935    }
 936
 937 }*/
 938
 939