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 */
  19
  20 #include "DataStructures.h"
  21 #include "AliRun.h"
  22 #include "AliDetector.h"
  23 #include "AliRICH.h"
  24 #include "AliRICHPoints.h"
  25 #include "AliRICHSegResV0.h"
  26 #include "AliRICHPatRec.h"
  27 #include "AliRICH.h"
  28 #include "AliRICHConst.h"
  29 #include "AliRICHPoints.h"
  30 #include "AliConst.h"
  31 #include "TParticle.h"
  32 #include "TMath.h"
  33 #include "TRandom.h"
  34 #include "TCanvas.h"
  35 #include "TH2.h"
  36
  37
  38 ClassImp(AliRICHPatRec)
  39 //___________________________________________
  40 AliRICHPatRec::AliRICHPatRec() : TObject()
  41 {
  42     //fChambers = 0;
  43 }
  44 //___________________________________________
  45 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
  46     : TObject()
  47 {
  48
  49 }
  50
  51 void AliRICHPatRec::PatRec()
  52 {
  53
  54   AliRICHChamber*       iChamber;
  55   AliRICHSegmentation*  segmentation;
  56
  57   Int_t ntracks, ndigits[7];
  58   Int_t itr, ich, i;
  59   Int_t GoodPhotons;
  60   Int_t x,y,q;
  61   Float_t rx,ry;
  62   Int_t nent,status;
  63
  64   Float_t gamma,MassCer,BetaCer;
  65
  66   Float_t rechit[5];
  67
  68   printf("PatRec started\n");
  69
  70   TCanvas *c1    = new TCanvas("c1","Alice RICH pad hits",50,10,700,700);
  71
  72   TH2F *ring     = new TH2F("ring","map",90,-30.,30.,90,-30.,30.);
  73   TH2F *ringband = new TH2F("ringband","map",90,-65.,65.,90,-65.,65.);
  74   TH1F *cerangle = new TH1F("cerangle","phot",300,0.45,0.75);
  75   TH1F *ceranglew= new TH1F("ceranglew","phot",300,0.45,0.75);
  76   TH1F *hough    = new TH1F("hough","hough",75,0.45,0.75);
  77   TH1F *mass     = new TH1F("mass","mass",100,50.,1050.);
  78
  79   AliRICH *RICH  = (AliRICH*)gAlice->GetDetector("RICH");
  80   TTree *TH = gAlice->TreeH();
  81
  82   ntracks =(Int_t) TH->GetEntries();
  83   //  ntracks = 1;
  84   for (itr=0; itr<ntracks; itr++) {
  85
  86     status = TrackParam(itr,ich);
  87     if(status==1) continue;
  88     printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
  89     //    ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
  90
  91     iChamber = &(RICH->Chamber(ich));
  92     segmentation=iChamber->GetSegmentationModel();
  93
  94     nent=(Int_t)gAlice->TreeD()->GetEntries();
  95     gAlice->TreeD()->GetEvent(nent-1);
  96     TClonesArray *Digits = RICH->DigitsAddress(ich);
  97     ndigits[ich] = Digits->GetEntriesFast();
  98     printf("ndigits %d in chamber %d\n",ndigits[ich],ich);
  99     AliRICHDigit *padI = 0;
 100
 101     GoodPhotons = 0;
 102
 103     for (Int_t dig=0;dig<ndigits[ich];dig++) {
 104       padI=(AliRICHDigit*) Digits->UncheckedAt(dig);
 105       x=padI->fPadX;
 106       y=padI->fPadY;
 107       q=padI->fSignal;
 108       segmentation->GetPadCxy(x,y,rx,ry);
 109
 110       fXpad = rx-fXshift;
 111       fYpad = ry-fYshift;
 112       fQpad = q;
 113
 114       ringband->Fill(x,y,1.);
 115       fCerenkovAnglePad = PhotonCerenkovAngle();
 116       if(fCerenkovAnglePad==-999) continue;
 117
 118       if(!PhotonInBand()) continue;
 119
 120       ring->Fill(fXpad,fYpad,1.);
 121       cerangle->Fill(fCerenkovAnglePad,1.);
 122
 123       GoodPhotons++;
 124       fEtaPhotons[GoodPhotons] = fCerenkovAnglePad;
 125     }
 126     fNumEtaPhotons = GoodPhotons;
 127
 128     BackgroundEstimation();
 129
 130     for(i=0;i<GoodPhotons;i++) {
 131       ceranglew->Fill(fEtaPhotons[i],fWeightPhotons[i]);
 132       //      printf(" Eta %f weight %f \n",fEtaPhotons[i],fWeightPhotons[i]);
 133     }
 134
 135     HoughResponse();
 136
 137     rechit[0] = 0;
 138     rechit[1]   = 0;
 139     rechit[2] = fThetaCerenkov;
 140     rechit[3] = 0;
 141     rechit[4] = 0;
 142
 143
 144     hough->Fill(fThetaCerenkov,1.);
 145
 146     RICH->AddRecHit(ich,rechit);
 147
 148     BetaCer = BetaCerenkov(1.29,fThetaCerenkov);
 149     gamma  = 1./sqrt(1.-pow(BetaCer,2));
 150     MassCer = fTrackMom/(BetaCer*gamma);
 151     //    printf(" mass %f \n",MassCer);
 152     mass->Fill(MassCer*1000,1.);
 153   }
 154
 155   gAlice->TreeR()->Fill();
 156   TClonesArray *fRec;
 157   for (i=0;i<7;i++) {
 158     fRec=RICH->RecHitsAddress(i);
 159     int ndig=fRec->GetEntriesFast();
 160     printf ("Chamber %d, rings %d\n",i,ndig);
 161   }
 162   RICH->ResetRecHits();
 163
 164
 165   c1->Divide(2,2);
 166   c1->cd(1);
 167   ring->Draw("box");
 168   c1->cd(2);
 169   ringband->Draw("box");
 170   c1->cd(3);
 171   cerangle->Draw();
 172   c1->cd(4);
 173   hough->Draw();
 174
 175 }
 176
 177
 178 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich)
 179 {
 180   // Get Local coordinates of track impact
 181
 182   AliRICHChamber*       iChamber;
 183   AliRICHSegmentation*  segmentation;
 184
 185   Float_t trackglob[3];
 186   Float_t trackloc[3];
 187   Float_t thetatr;
 188   Float_t phitr;
 189   Float_t iloss;
 190   Float_t part;
 191   Float_t pX, pY, pZ;
 192
 193   printf("Calling TrackParam\n");
 194
 195     gAlice->ResetHits();
 196     TTree *TH = gAlice->TreeH();
 197     TH->GetEvent(itr);
 198
 199     AliRICH *RICH  = (AliRICH*)gAlice->GetDetector("RICH");
 200     AliRICHHit* mHit=(AliRICHHit*)RICH->FirstHit(-1);
 201     if(mHit==0) return 1;
 202     ich = mHit->fChamber-1;
 203     trackglob[0] = mHit->fX;
 204     trackglob[1] = mHit->fY;
 205     trackglob[2] = mHit->fZ;
 206     pX = mHit->fMomX;
 207     pY = mHit->fMomY;
 208     pZ = mHit->fMomZ;
 209     fTrackMom = sqrt(pow(pX,2)+pow(pY,2)+pow(pZ,2));
 210     thetatr = (180 - mHit->fTheta)*(Float_t)kDegrad;
 211     phitr = mHit->fPhi*(Float_t)kDegrad;
 212     iloss = mHit->fLoss;
 213     part  = mHit->fParticle;
 214
 215     iChamber = &(RICH->Chamber(ich));
 216     iChamber->GlobaltoLocal(trackglob,trackloc);
 217
 218     segmentation=iChamber->GetSegmentationModel();
 219
 220     // retrieve geometrical params
 221
 222     AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();
 223
 224     fRw   = fGeometry->GetFreonThickness();
 225     fQw   = fGeometry->GetQuartzThickness();
 226     fTgap = fGeometry->GetGapThickness()
 227             + fGeometry->GetProximityGapThickness();
 228
 229     Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
 230     fTrackLoc[0] = apar*cos(phitr);
 231     fTrackLoc[1] = apar*sin(phitr);
 232     fTrackLoc[2] = fRw + fQw + fTgap;
 233     fTrackTheta  = thetatr;
 234     fTrackPhi    = phitr;
 235
 236     fXshift = trackloc[0] - fTrackLoc[0];
 237     fYshift = trackloc[2] - fTrackLoc[1];
 238
 239     return 0;
 240 }
 241
 242 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
 243                                                  Float_t phiphot)
 244 {
 245   Float_t nquartz = 1.585;
 246   Float_t ngas    = 1.;
 247   Float_t nfreon  = 1.295;
 248   Float_t value;
 249
 250   //  printf("Calling EstimationLimits\n");
 251
 252   Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
 253   Float_t b1 = (fRw-fEmissPoint)*tan(lim);
 254   Float_t b2 = fQw / sqrt(pow(nquartz,2)-pow(nfreon*sin(lim),2));
 255   Float_t b3 = fTgap / sqrt(pow(ngas,2)-pow(nfreon*sin(lim),2));
 256   Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
 257   value = pow(radius,2)
 258     -pow((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
 259     -pow((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
 260   return value;
 261 }
 262
 263
 264 Float_t AliRICHPatRec::PhotonCerenkovAngle()
 265 {
 266   // Cherenkov pad angle reconstruction
 267
 268   Float_t radius;
 269   Float_t cherMin = 0;
 270   Float_t cherMax = 0.8;
 271   Float_t phiphot;
 272   Float_t eps = 0.0001;
 273   Int_t   niterEmiss = 0;
 274   Int_t   niterEmissMax = 0;
 275   Float_t x1,x2,x3,p1,p2,p3;
 276   Float_t argY,argX;
 277   Int_t niterFun;
 278
 279   //  printf("Calling PhotonCerenkovAngle\n");
 280
 281   radius = sqrt(pow(fTrackLoc[0]-fXpad,2)+pow(fTrackLoc[1]-fYpad,2));
 282   fEmissPoint = fRw/2.;  //Start value of EmissionPoint
 283
 284   while(niterEmiss<=niterEmissMax) {
 285
 286     niterFun = 0;
 287     argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
 288     argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
 289     phiphot = atan2(argY,argX);
 290     p1 = EstimationAtLimits(cherMin,radius,phiphot);
 291     p2 = EstimationAtLimits(cherMax,radius,phiphot);
 292     if(p1*p2>0)
 293       {
 294         //        printf("PhotonCerenkovAngle failed\n");
 295         return -999;
 296       }
 297
 298     //start to find the Cherenkov pad angle
 299     x1 = cherMin;
 300     x2 = cherMax;
 301     x3 = (x1+x2)/2.;
 302     p3 = EstimationAtLimits(x3,radius,phiphot);
 303     while(TMath::Abs(p3)>eps){
 304       if(p1*p3<0) x2 = x3;
 305       if(p1*p3>0) {
 306         x1 = x3;
 307         p1 = EstimationAtLimits(x1,radius,phiphot);
 308       }
 309       x3 = (x1+x2)/2.;
 310       p3 = EstimationAtLimits(x3,radius,phiphot);
 311       niterFun++;
 312
 313       if(niterFun>=1000) {
 314         //      printf(" max iterations in PhotonCerenkovAngle\n");
 315         return x3;
 316       }
 317     }
 318     //    printf("niterFun %i \n",niterFun);
 319     niterEmiss++;
 320     if (niterEmiss != niterEmissMax+1) EmissionPoint();
 321   }
 322   /*
 323    printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
 324                          phiphot,fXpad,fYpad,fEmissPoint);
 325   */
 326
 327   return x3;
 328
 329 }
 330
 331
 332 void AliRICHPatRec::EmissionPoint()
 333 {
 334   Float_t AbsorLength=7.83*fRw; //absorption length in the freon (cm)
 335   // 7.83 = -1/ln(T0) where
 336   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
 337   Float_t PhotonLength, PhotonLengthMin, PhotonLengthMax;
 338
 339   PhotonLength=exp(-fRw/(AbsorLength*cos(fCerenkovAnglePad)));
 340   PhotonLengthMin=fRw*PhotonLength/(1.-PhotonLength);
 341   PhotonLengthMax=AbsorLength*cos(fCerenkovAnglePad);
 342   fEmissPoint = fRw + PhotonLengthMin - PhotonLengthMax;
 343
 344 }
 345
 346 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
 347 {
 348   printf("Calling PhotonSelection\n");
 349 }
 350
 351 void AliRICHPatRec::BackgroundEstimation()
 352 {
 353   Float_t StepEta   = 0.001;
 354   Float_t EtaMinBkg = 0.72;
 355   Float_t EtaMaxBkg = 0.75;
 356   Float_t EtaMin    = 0.;
 357   Float_t EtaMax    = 0.75;
 358   Float_t ngas      = 1.;
 359   Float_t nfreon    = 1.295;
 360
 361   Float_t EtaStepMin,EtaStepMax,EtaStepAvg;
 362   Int_t i,ip,nstep;
 363   Int_t NumPhotBkg, NumPhotonStep;
 364   Float_t FunBkg,AreaBkg,NormBkg;
 365   Float_t DensityBkg,StoreBkg,NumStore;
 366   Float_t ThetaSig;
 367
 368   NumPhotBkg = 0;
 369   AreaBkg = 0.;
 370
 371   nstep = (int)((EtaMaxBkg-EtaMinBkg)/StepEta);
 372
 373   for (i=0;i<fNumEtaPhotons;i++) {
 374
 375     if(fEtaPhotons[i]>EtaMinBkg && fEtaPhotons[i]<EtaMaxBkg) {
 376       NumPhotBkg++;
 377     }
 378   }
 379   if (NumPhotBkg == 0) {
 380      for (i=0;i<fNumEtaPhotons;i++) {
 381         fWeightPhotons[i] = 1.;
 382      }
 383     return;
 384   }
 385
 386   //  printf(" NumPhotBkg %i ",NumPhotBkg);
 387
 388   for (i=0;i<nstep;i++) {
 389     EtaStepMin = EtaMinBkg + (Float_t)(i)*StepEta;
 390     EtaStepMax = EtaMinBkg + (Float_t)(i+1)*StepEta;
 391     EtaStepAvg = 0.5*(EtaStepMax + EtaStepMin);
 392     /*
 393     FunBkg = tan(EtaStepAvg)*pow((1.+pow(tan(EtaStepAvg),2)),
 394                                   5.52)-7.803 + 22.02*tan(EtaStepAvg);
 395     */
 396     ThetaSig = asin(nfreon/ngas*sin(EtaStepAvg));
 397     FunBkg = tan(ThetaSig)*(1.+pow(tan(ThetaSig),2))*nfreon
 398        /ngas*cos(EtaStepAvg)/cos(ThetaSig);
 399     AreaBkg += StepEta*FunBkg;
 400   }
 401
 402   DensityBkg = 0.95*(Float_t)(NumPhotBkg)/AreaBkg;
 403   //  printf(" DensityBkg %f \n",DensityBkg);
 404
 405   nstep = (int)((EtaMax-EtaMin)/StepEta);
 406   StoreBkg = 0.;
 407   NumStore = 0;
 408   for (i=0;i<nstep;i++) {
 409     EtaStepMin = EtaMinBkg + (Float_t)(i)*StepEta;
 410     EtaStepMax = EtaMinBkg + (Float_t)(i+1)*StepEta;
 411     EtaStepAvg = 0.5*(EtaStepMax + EtaStepMin);
 412     /*
 413     FunBkg = tan(EtaStepAvg)*pow((1.+pow(tan(EtaStepAvg),2)),
 414                                   5.52)-7.803 + 22.02*tan(EtaStepAvg);
 415     */
 416
 417     ThetaSig = asin(nfreon/ngas*sin(EtaStepAvg));
 418     FunBkg = tan(ThetaSig)*(1.+pow(tan(ThetaSig),2))*nfreon
 419        /ngas*cos(EtaStepAvg)/cos(ThetaSig);
 420
 421     AreaBkg = StepEta*FunBkg;
 422     NormBkg = DensityBkg*AreaBkg;
 423     NumPhotonStep = 0;
 424     for (ip=0;ip<fNumEtaPhotons;ip++) {
 425       if(fEtaPhotons[ip]>EtaStepMin && fEtaPhotons[ip]<EtaStepMax) {
 426         NumPhotonStep++;
 427       }
 428     }
 429     if (NumPhotonStep == 0) {
 430       StoreBkg += NormBkg;
 431       NumStore++;
 432       if (NumStore>50) {
 433         NumStore = 0;
 434         StoreBkg = 0.;
 435       }
 436     }
 437     if (NumPhotonStep == 0) continue;
 438     for (ip=0;ip<fNumEtaPhotons;ip++) {
 439       if(fEtaPhotons[ip]>EtaStepMin && fEtaPhotons[ip]<EtaStepMax) {
 440         NormBkg +=StoreBkg;
 441         StoreBkg = 0;
 442         NumStore = 0;
 443         fWeightPhotons[ip] = 1. - NormBkg/(Float_t)(NumPhotonStep);
 444         /*
 445         printf(" NormBkg %f NumPhotonStep %i fW %f \n",
 446                NormBkg, NumPhotonStep, fWeightPhotons[ip]);
 447         */
 448         if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
 449       }
 450     }
 451   }
 452 }
 453
 454
 455 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
 456 {
 457   printf("Calling FlagPhotons\n");
 458 }
 459
 460
 461 //////////////////////////////////////////
 462
 463
 464
 465
 466
 467 Int_t AliRICHPatRec::PhotonInBand()
 468 {
 469   //0=label for parameters giving internal band ellipse
 470   //1=label for parameters giving external band ellipse
 471
 472   Float_t imp[2], mass[2], Energ[2], beta[2];
 473   Float_t EmissPointLength[2];
 474   Float_t E1, E2, F1, F2;
 475   Float_t nfreon[2], nquartz[2];
 476   Int_t times;
 477
 478
 479   Float_t phpad, thetacer[2];
 480   Float_t bandradius[2], padradius;
 481
 482   imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
 483   imp[1] = 1.2;
 484
 485   mass[0] = 0.938; //proton mass
 486   mass[1] = 0.139; //pion mass
 487
 488   EmissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
 489   EmissPointLength[1] = 0.;//at the end of radiator
 490
 491   //parameters to calculate freon window refractive index vs. energy
 492   Float_t a = 1.177;
 493   Float_t b = 0.0172;
 494
 495   //parameters to calculate quartz window refractive index vs. energy
 496   /*
 497   Energ[0]  = 5.6;
 498   Energ[1]  = 7.7;
 499   */
 500   Energ[0]  = 5.0;
 501   Energ[1]  = 8.0;
 502   E1 = 10.666;
 503   E2  = 18.125;
 504   F1  = 46.411;
 505   F2  = 228.71;
 506
 507
 508   phpad = PhiPad();
 509
 510   for (times=0; times<=1; times++) {
 511
 512     nfreon[times]   = a+b*Energ[times];
 513
 514     nquartz[times] = sqrt(1+(F1/(pow(E1,2)-pow(Energ[times],2)))+
 515                           (F2/(pow(E2,2)-pow(Energ[times],2))));
 516
 517     beta[times]  = imp[times]/sqrt(pow(imp[times],2)+pow(mass[times],2));
 518
 519     thetacer[times] =  CherenkovAngle( nfreon[times], beta[times]);
 520
 521     bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
 522                                         EmissPointLength[times],
 523                                         thetacer[times], phpad);
 524   }
 525
 526   bandradius[0] -= 1.6;
 527   bandradius[1] += 1.6;
 528   padradius = sqrt(pow(fXpad,2)+pow(fYpad,2));
 529   //  printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
 530
 531   if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
 532   return 0;
 533 }
 534
 535 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
 536                        Float_t EmissPointLength, Float_t thetacer,
 537                        Float_t phpad)
 538 {
 539   Float_t DistanceValue;
 540
 541   TVector3 RadExitPhot(1,1,1);//photon impact at the radiator exit with respect
 542   //to local reference sistem with the origin in the MIP entrance
 543
 544   TVector3 VectEmissPointLength(1,1,1);
 545   Float_t MagEmissPointLenght;
 546
 547   TVector3 RadExitPhot2(1,1,1);//photon impact at the radiator exit with respect
 548   Float_t MagRadExitPhot2;
 549   //to a reference sistem with origin in the photon emission point and
 550   //axes parallel to the MIP reference sistem
 551
 552   TVector3 QuarExitPhot(1,1,1);//photon impact at the quartz exit with respect
 553   Float_t MagQuarExitPhot;
 554   //
 555   TVector3 GapExitPhot(1,1,1) ;
 556   Float_t MagGapExitPhot;
 557   //
 558   TVector3 fPhotocatExitPhot(1,1,1);
 559   Double_t theta2;
 560   Double_t thetarad , phirad ;
 561   Double_t thetaquar, phiquar;
 562   Double_t thetagap , phigap ;
 563
 564   Float_t ngas    = 1.;
 565
 566   MagEmissPointLenght =  EmissPointLength/cos(fTrackTheta);
 567
 568   VectEmissPointLength.SetMag(MagEmissPointLenght);
 569   VectEmissPointLength.SetTheta(fTrackTheta);
 570   VectEmissPointLength.SetPhi(fTrackPhi);
 571
 572
 573   RadExitPhot2.SetTheta(thetacer);
 574   RadExitPhot2.SetPhi(phpad);
 575
 576
 577   TRotation r1;
 578   TRotation r2;
 579   TRotation r;
 580
 581   r1. RotateY(fTrackTheta);
 582   r2. RotateZ(fTrackPhi);
 583
 584
 585
 586   r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
 587   //following by a rotation about the z axis by MIP phi incidence angle;
 588
 589
 590   RadExitPhot2    = r * RadExitPhot2;
 591   theta2          = RadExitPhot2.Theta();
 592   MagRadExitPhot2 = (fRw -  VectEmissPointLength(2))/cos(theta2);
 593   RadExitPhot2.SetMag(MagRadExitPhot2);
 594
 595
 596   RadExitPhot = VectEmissPointLength + RadExitPhot2;
 597   thetarad    = RadExitPhot.Theta();
 598
 599   phirad  =  RadExitPhot.Phi(); //check on the original file //
 600
 601   thetaquar   = SnellAngle( nfreon, nquartz, theta2);
 602   phiquar     = RadExitPhot2.Phi();
 603   if(thetaquar == 999.) return thetaquar;
 604   MagQuarExitPhot    = fQw/cos(thetaquar);
 605   QuarExitPhot.SetMag( MagQuarExitPhot);
 606   QuarExitPhot.SetTheta(thetaquar);
 607   QuarExitPhot.SetPhi(phiquar);
 608
 609   thetagap = SnellAngle( nquartz, ngas, thetaquar);
 610   phigap   = phiquar;
 611   if(thetagap == 999.) return thetagap;
 612   MagGapExitPhot    = fTgap/cos(thetagap);
 613   GapExitPhot.SetMag( MagGapExitPhot);
 614   GapExitPhot.SetTheta(thetagap);
 615   GapExitPhot.SetPhi(phigap);
 616
 617   fPhotocatExitPhot =  RadExitPhot + QuarExitPhot + GapExitPhot;
 618
 619   DistanceValue = sqrt(pow(fPhotocatExitPhot(0),2)
 620                            +pow(fPhotocatExitPhot(1),2));
 621   return  DistanceValue ;
 622 }
 623
 624 Float_t AliRICHPatRec::PhiPad()
 625 {
 626   Float_t zpad;
 627   Float_t thetapad, phipad;
 628   Float_t thetarot, phirot;
 629
 630   zpad = fRw + fQw + fTgap;
 631
 632   TVector3 PhotonPad(fXpad, fYpad, zpad);
 633   thetapad = PhotonPad.Theta();
 634   phipad = PhotonPad.Phi();
 635
 636   TRotation r1;
 637   TRotation r2;
 638   TRotation r;
 639
 640   thetarot = - fTrackTheta;
 641   phirot   = - fTrackPhi;
 642   r1. RotateZ(phirot);
 643   r2. RotateY(thetarot);
 644
 645   r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
 646   //following by a rotation about the y axis by MIP -theta incidence angle;
 647
 648   PhotonPad  = r * PhotonPad;
 649
 650   phipad = PhotonPad.Phi();
 651
 652   return phipad;
 653 }
 654
 655 Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
 656 {
 657   Float_t sinrefractangle;
 658   Float_t refractangle;
 659
 660   sinrefractangle = (n1/n2)*sin(theta1);
 661
 662   if(sinrefractangle>1.) {
 663      refractangle = 999.;
 664      return refractangle;
 665   }
 666
 667   refractangle = asin(sinrefractangle);
 668   return refractangle;
 669 }
 670
 671 Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
 672 {
 673   Float_t thetacer;
 674
 675   if((n*beta)<1.) {
 676     thetacer = 999.;
 677     return thetacer;
 678   }
 679
 680   thetacer = acos (1./(n*beta));
 681   return thetacer;
 682 }
 683
 684 Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
 685 {
 686   Float_t beta;
 687
 688   beta = 1./(n*cos(theta));
 689   return beta;
 690 }
 691
 692
 693
 694
 695 void AliRICHPatRec::HoughResponse()
 696
 697 {
 698   int           bin=0;
 699   int           bin1=0;
 700   int           bin2=0;
 701   int           i, j, k, NcorrBand;
 702   int           EtaBin = 750;
 703   float         HCS[750];
 704   float         angle, ThetaCerMean;
 705
 706   float         EtaPeak[30];
 707   float         EtaMin = 0.00;
 708   float         EtaMax = 0.75;
 709   float         StepEta = 0.001;
 710   float         WindowEta = 0.040;
 711
 712   int           Nbin;
 713
 714   float EtaPeakPos  = -1;
 715   Int_t   EtaPeakCount = -1;
 716
 717   ThetaCerMean   = 0.;
 718   fThetaCerenkov = 0.;
 719
 720   Nbin = (int)(0.5+EtaMax/(StepEta));
 721   NcorrBand = (int)(0.5+ WindowEta/(2 * StepEta));
 722   memset ((void *)HCS, 0, EtaBin*sizeof(int));
 723
 724   for (k=0; k< fNumEtaPhotons; k++) {
 725
 726     angle = fEtaPhotons[k];
 727
 728     if (angle>=EtaMin && angle<= EtaMax) {
 729       bin = (int)(0.5+angle/(StepEta));
 730       bin1= bin-NcorrBand;
 731       bin2= bin+NcorrBand;
 732       if (bin1<0)    bin1=0;
 733       if (bin2>Nbin) bin2=Nbin;
 734
 735       for (j=bin1; j<bin2; j++) {
 736         HCS[j] += fWeightPhotons[k];
 737       }
 738
 739       ThetaCerMean += angle;
 740     }
 741   }
 742
 743  ThetaCerMean /= fNumEtaPhotons;
 744
 745   HoughFiltering(HCS);
 746
 747   for (bin=0; bin <Nbin; bin++) {
 748     angle = (bin+0.5) * (StepEta);
 749     if (HCS[bin] && HCS[bin] > EtaPeakPos) {
 750       EtaPeakCount = 0;
 751       EtaPeakPos = HCS[bin];
 752       EtaPeak[0]=angle;
 753     }
 754     else {
 755       if (HCS[bin] == EtaPeakPos) {
 756         EtaPeak[++EtaPeakCount] = angle;
 757       }
 758     }
 759   }
 760
 761   for (i=0; i<EtaPeakCount+1; i++) {
 762     fThetaCerenkov += EtaPeak[i];
 763   }
 764   if (EtaPeakCount>=0) {
 765     fThetaCerenkov /= EtaPeakCount+1;
 766     fThetaPeakPos = EtaPeakPos;
 767   }
 768 }
 769
 770
 771 void AliRICHPatRec::HoughFiltering(float HCS[])
 772 {
 773    float HCS_filt[750];
 774    float K[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
 775    int nx, i, nx_dx;
 776    int sizeHCS;
 777    int Nbin;
 778
 779    int   EtaBin = 750;
 780    float EtaMax = 0.75;
 781    float StepEta = 0.001;
 782
 783    Nbin =  (int)(1+EtaMax/StepEta);
 784    sizeHCS = EtaBin*sizeof(float);
 785
 786    memset ((void *)HCS_filt, 0, sizeHCS);
 787
 788    for (nx = 0; nx < Nbin; nx++) {
 789       for (i = 0; i < 5; i++)   {
 790         nx_dx = nx + (i-2);
 791         if (nx_dx> -1 && nx_dx<Nbin)
 792              HCS_filt[nx] +=  HCS[nx_dx] * K[i];
 793       }
 794    }
 795
 796    for (nx = 0; nx < Nbin; nx++) {
 797      HCS[nx] = HCS_filt[nx];
 798    }
 799 }
 800
 801 Float_t AliRICHPatRec::CherenkovRingDrawing(Float_t fixedthetacer)
 802
 803 {
 804
 805 //to draw Cherenkov ring by known Cherenkov angle
 806
 807    Int_t nmaxdegrees, nstepdegrees;
 808    Float_t phpad, thetacer;
 809    Float_t nfreonave, nquartzave;
 810    Float_t AveEnerg;
 811    Float_t Energ[2];
 812    Float_t E1, E2, F1, F2;
 813    Float_t bandradius;
 814    Float_t CoordPadRing;
 815
 816 //parameters to calculate freon window refractive index vs. energy
 817    Float_t a = 1.177;
 818    Float_t b = 0.0172;
 819
 820 //parameters to calculate quartz window refractive index vs. energy
 821 /*
 822    Energ[0]  = 5.6;
 823    Energ[1]  = 7.7;
 824 */
 825    Energ[0]  = 5.0;
 826    Energ[1]  = 8.0;
 827    E1  = 10.666;
 828    E2  = 18.125;
 829    F1  = 46.411;
 830    F2  = 228.71;
 831
 832
 833    nmaxdegrees = 360;
 834
 835    nstepdegrees = 36;
 836
 837    for (phpad=0; phpad<nmaxdegrees;phpad++) {
 838
 839      AveEnerg =  (Energ[0]+Energ[1])/2.;
 840
 841      nfreonave  = a+b*AveEnerg;
 842      nquartzave = sqrt(1+(F1/(pow(E1,2)-pow(AveEnerg,2)))+
 843                          (F2/(pow(E2,2)-pow(AveEnerg,2))));
 844
 845      thetacer =  fixedthetacer;
 846
 847      bandradius = DistanceFromMip(nfreonave, nquartzave,
 848                                    fEmissPoint,thetacer, phpad);
 849
 850      CoordPadRing=fPhotocatExitPhot;
 851
 852      phpad = (nmaxdegrees/nstepdegrees)*phpad;
 853
 854      return CoordPadRing;
 855                                                                                     }
 856  }