RICH/AliRICHRecon.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 //                                                                      //
  18 // AliRICHRecon                                                         //
  19 //                                                                      //
  20 // RICH class to perfom pattern recognition based on Hough transfrom    //
  21 //                                                                      //
  22 //////////////////////////////////////////////////////////////////////////
  23
  24
  25 #include "AliRICH.h"
  26 #include "AliRICHRecon.h"
  27 #include "AliRICHParam.h"
  28 #include "AliRICHChamber.h"
  29 #include <AliLoader.h>
  30 #include <AliRun.h>
  31 #include <AliStack.h>
  32 #include <Riostream.h>
  33 #include <TParticle.h>
  34 #include <TH2.h>
  35 #include <TMath.h>
  36 #include <TRandom.h>
  37 #include <TMinuit.h>
  38 #include <TNtuple.h>
  39 #include <TMath.h>
  40 #include <TRotation.h>
  41 #include <TVector3.h>
  42 #include <TCanvas.h>
  43
  44 #define NPointsOfRing 201
  45
  46 TMinuit *gAliRICHminuit ;
  47
  48 void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  49 //__________________________________________________________________________________________________
  50 AliRICHRecon::AliRICHRecon(const char*name, const char*title)
  51              :TTask(name,title)
  52 {
  53   // main ctor
  54   fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
  55   fDTheta       = 0.001;   fWindowWidth  = 0.060;
  56   fNpadX         = AliRICHParam::NpadsY();
  57   fNpadY         = AliRICHParam::NpadsX();
  58   fPadSizeX      = AliRICHParam::PadSizeY();
  59   fPadSizeY      = AliRICHParam::PadSizeX();
  60   fRadiatorWidth = AliRICHParam::FreonThickness();
  61   fQuartzWidth   = AliRICHParam::QuartzThickness();
  62   fGapWidth      = AliRICHParam::RadiatorToPads();
  63   fXmin         = -AliRICHParam::PcSizeY()/2.;
  64   fXmax         =  AliRICHParam::PcSizeY()/2.;
  65   fYmin         = -AliRICHParam::PcSizeX()/2.;
  66   fYmax         =  AliRICHParam::PcSizeX()/2.;
  67   fRich = (AliRICH*)gAlice->GetDetector("RICH");
  68   fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos");
  69   if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750);
  70   fNtuple=new TNtuple("hn","ntuple",
  71 "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ");
  72 }
  73 //__________________________________________________________________________________________________
  74 void AliRICHRecon::StartProcessEvent()
  75 {
  76   //start to process for pattern recognition
  77
  78   Float_t trackThetaStored    = 0;
  79   Float_t trackPhiStored      = 0;
  80   Float_t thetaCerenkovStored = 0;
  81   Int_t houghPhotonsStored    = 0;
  82
  83   SetFreonScaleFactor(0.994);
  84
  85   if(fIsDISPLAY)
  86     {
  87       DrawEvent(0);
  88 //      Waiting();
  89     }
  90
  91     Rich()->GetLoader()->LoadHits();
  92     Rich()->GetLoader()->LoadRecPoints();
  93     Rich()->GetLoader()->LoadDigits();
  94     gAlice->GetRunLoader()->LoadHeader();
  95     gAlice->GetRunLoader()->LoadKinematics();
  96
  97     Rich()->GetLoader()->TreeR()->GetEntry(0);
  98
  99     Float_t clusX[7][500],clusY[7][500];
 100     Int_t clusQ[7][500],clusMul[7][500];
 101     Int_t nClusters[7];
 102
 103     for (Int_t ich=0;ich<7;ich++) {
 104       nClusters[ich] = Rich()->Clusters(ich+1)->GetEntries();
 105       for(Int_t k=0;k<nClusters[ich];k++) {
 106         AliRICHcluster *pCluster = (AliRICHcluster *)Rich()->Clusters(ich+1)->At(k);
 107         clusX[ich][k] = pCluster->X();
 108         clusY[ich][k] = pCluster->Y();
 109         clusQ[ich][k] = pCluster->Q();
 110         clusMul[ich][k] = pCluster->Size();
 111 //        pCluster->Print();
 112       }
 113     }
 114
 115     Int_t nPrimaries = (Int_t)Rich()->GetLoader()->TreeH()->GetEntries();
 116
 117     cout << " N. primaries " << nPrimaries << endl;
 118
 119     for(Int_t i=0;i<nPrimaries;i++){
 120
 121       Rich()->GetLoader()->TreeH()->GetEntry(i);
 122
 123 //      Rich()->Hits()->Print();
 124       Int_t iPrim = 0;
 125
 126       AliRICHhit* pHit=0;
 127
 128       for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) {
 129
 130         pHit = (AliRICHhit*)Rich()->Hits()->At(j);
 131         if(pHit->GetTrack() < nPrimaries) break;
 132         iPrim++;
 133       }
 134
 135       cout << " iPrim " << iPrim << " pHit " << pHit << endl;
 136
 137       if (!pHit) return;
 138
 139 //      pHit->Print();
 140
 141       TParticle *pParticle = gAlice->GetRunLoader()->Stack()->Particle(pHit->GetTrack());
 142       Float_t pmod     = pParticle->P();
 143       Float_t pt       = pParticle->Pt();
 144       Float_t trackEta = pParticle->Eta();
 145       Int_t q          = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge());
 146
 147 //      pParticle->Print();
 148
 149       cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl;
 150
 151       SetTrackMomentum(pmod);
 152       SetTrackPt(pt);
 153       SetTrackEta(trackEta);
 154       SetTrackCharge(q);
 155
 156       TVector3 pLocal(0,0,0);//?????
 157
 158       TVector2 primLocal =Rich()->C(pHit->C())->Glob2Loc(pHit->InX3());
 159
 160 //      Float_t pmodFreo = pLocal.Mag();
 161       Float_t trackTheta = pLocal.Theta();
 162       Float_t trackPhi = pLocal.Phi();
 163
 164 //      cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl;
 165
 166       SetTrackTheta(trackTheta);
 167       SetTrackPhi(trackPhi);
 168
 169       Int_t maxInd = 0;
 170       Float_t minDist =  999.;
 171
 172 //      cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl;
 173
 174       for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++)
 175         {
 176           Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j];
 177           Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j];
 178
 179
 180           Float_t diff = sqrt(diffx*diffx + diffy*diffy);
 181
 182           if(diff < minDist)
 183             {
 184               minDist = diff;
 185               maxInd = j;
 186             }
 187
 188         }
 189
 190       Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd];
 191       Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd];
 192
 193       cout << " diffx " << diffx << " diffy " << diffy << endl;
 194
 195
 196       SetMipIndex(maxInd);
 197       SetTrackIndex(i);
 198
 199       Float_t shiftX = 0;//primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X(); ?????
 200       Float_t shiftY = 0;//primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y(); ?????
 201
 202       SetShiftX(shiftX);
 203       SetShiftY(shiftY);
 204
 205       Float_t *pclusX = &clusX[pHit->Chamber()-1][0];
 206       Float_t *pclusY = &clusY[pHit->Chamber()-1][0];
 207
 208       SetCandidatePhotonX(pclusX);
 209       SetCandidatePhotonY(pclusY);
 210       SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]);
 211
 212       Int_t qch = clusQ[pHit->Chamber()-1][maxInd];
 213
 214
 215       if(minDist < 3.0 && qch > 120 && maxInd !=0)
 216         {
 217
 218           if(fIsBACKGROUND)
 219             {
 220
 221               Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964);
 222               Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964);
 223               SetShiftX(xrndm);
 224               SetShiftY(yrndm);
 225
 226             }
 227
 228           PatRec();
 229
 230           trackThetaStored = GetTrackTheta();
 231           trackPhiStored = GetTrackPhi();
 232           thetaCerenkovStored = GetThetaCerenkov();
 233           houghPhotonsStored = GetHoughPhotons();
 234
 235           Int_t diffNPhotons = 999;
 236           Int_t nsteps = 0;
 237           Float_t diffTrackTheta = 999.;
 238           Float_t diffTrackPhi   = 999.;
 239
 240           while(fIsMINIMIZER && GetHoughPhotons() > 2
 241                             && diffNPhotons !=0
 242                             && diffTrackTheta > 0.0001
 243                             && nsteps < 10)
 244             {
 245
 246               Int_t   houghPhotonsBefore  = GetHoughPhotons();
 247
 248               Float_t trackThetaBefore = GetTrackTheta();
 249               Float_t trackPhiBefore   = GetTrackPhi();
 250
 251               Minimization();
 252
 253               PatRec();
 254
 255               diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons());
 256
 257               Float_t trackThetaAfter = GetTrackTheta();
 258               Float_t trackPhiAfter   = GetTrackPhi();
 259
 260               diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore);
 261               diffTrackPhi   = TMath::Abs(trackPhiAfter - trackPhiBefore);
 262
 263               if(fDebug)
 264               cout << " houghPhotonsBefore " << houghPhotonsBefore
 265                    << " GetHoughPhotons()  " << GetHoughPhotons();
 266
 267               nsteps++;
 268             }
 269
 270           SetFittedThetaCerenkov(GetThetaCerenkov());
 271           SetFittedHoughPhotons(GetHoughPhotons());
 272
 273           SetTrackTheta(trackThetaStored);
 274           SetTrackPhi(trackPhiStored);
 275           SetThetaCerenkov(thetaCerenkovStored);
 276           SetHoughPhotons(houghPhotonsStored);
 277
 278           SetMinDist(minDist);
 279
 280           FillHistograms();
 281
 282           if(fIsDISPLAY) DrawEvent(1);
 283
 284           Waiting();
 285
 286         }
 287     }
 288   if(fIsDISPLAY) fDisplay->Print("display.ps");
 289 }//StartProcessEvent()
 290 //__________________________________________________________________________________________________
 291 void AliRICHRecon::EndProcessEvent()
 292 {
 293   // function called at the end of the event loop
 294
 295   fOutFile->Write();
 296   fOutFile->Close();
 297 }
 298 //__________________________________________________________________________________________________
 299 void AliRICHRecon::PatRec()
 300 {
 301   //pattern recognition method based on Hough transform
 302
 303
 304   Float_t trackTheta = GetTrackTheta();
 305   Float_t trackPhi   = GetTrackPhi();
 306   Float_t pmod       = GetTrackMomentum();
 307   Int_t iMipIndex   = GetMipIndex();
 308
 309   Bool_t kPatRec = kFALSE;
 310
 311   Int_t candidatePhotons = 0;
 312
 313   Float_t shiftX = GetShiftX();
 314   Float_t shiftY = GetShiftY();
 315
 316   Float_t* candidatePhotonX = GetCandidatePhotonX();
 317   Float_t* candidatePhotonY = GetCandidatePhotonY();
 318
 319   Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
 320
 321   if(fDebug) cout << " n " << candidatePhotonsNumber << endl;
 322
 323   SetThetaCerenkov(999.);
 324   SetHoughPhotons(0);
 325   SetHoughPhotonsNorm(0);
 326   SetHoughRMS(999.);
 327
 328   for (Int_t j=0; j < candidatePhotonsNumber; j++)
 329     {
 330
 331       SetPhotonIndex(j);
 332
 333       SetPhotonFlag(0);
 334       SetPhotonEta(-999.);
 335       SetPhotonWeight(0.);
 336
 337       if (j == iMipIndex) continue;
 338
 339
 340       if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */
 341
 342       Float_t xtoentr = candidatePhotonX[j] - shiftX;
 343       Float_t ytoentr = candidatePhotonY[j] - shiftY;
 344
 345       //      Float_t chargehit = fHits_charge[j];
 346       //      if(chargehit > 150) continue;
 347
 348       SetEntranceX(xtoentr);
 349       SetEntranceY(ytoentr);
 350
 351       FindPhiPoint();
 352
 353       Int_t photonStatus = PhotonInBand();
 354
 355       if(fDebug)
 356          {
 357             cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl;
 358             cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl;
 359          }
 360
 361       if(photonStatus == 0) continue;
 362
 363       SetPhotonFlag(1);
 364
 365       FindThetaPhotonCerenkov();
 366
 367       Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
 368
 369       if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl;
 370
 371       SetPhotonEta(thetaPhotonCerenkov);
 372
 373       candidatePhotons++;
 374
 375
 376     }
 377
 378   if(candidatePhotons >= 1) kPatRec = kTRUE;
 379
 380   if(!kPatRec) return;
 381     {
 382        SetThetaCerenkov(999.);
 383        SetHoughPhotons(0);
 384     }
 385   SetPhotonsNumber(candidatePhotonsNumber);
 386
 387   HoughResponse();
 388
 389   fNrings++;
 390
 391   FlagPhotons();
 392   Int_t nPhotonHough = GetHoughPhotons();
 393
 394   if(nPhotonHough < 1)
 395     {
 396       SetThetaCerenkov(999.);
 397       SetHoughPhotonsNorm(0.);
 398       return;
 399     }
 400
 401   if(fIsWEIGHT) FindWeightThetaCerenkov();
 402
 403   Float_t thetaCerenkov = GetThetaCerenkov();
 404
 405   SetThetaOfRing(thetaCerenkov);
 406   FindAreaAndPortionOfRing();
 407
 408   Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
 409   SetHoughPhotonsNorm(nPhotonHoughNorm);
 410
 411   // Calculate the area where the photon are accepted...
 412
 413   Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
 414   SetThetaOfRing(thetaInternal);
 415   FindAreaAndPortionOfRing();
 416   Float_t internalArea = GetAreaOfRing();
 417
 418   Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
 419   SetThetaOfRing(thetaExternal);
 420   FindAreaAndPortionOfRing();
 421   Float_t externalArea = GetAreaOfRing();
 422
 423   Float_t houghArea = externalArea - internalArea;
 424
 425   SetHoughArea(houghArea);
 426
 427   if(fDebug)
 428     {
 429       cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl;
 430       cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl;
 431
 432
 433       cout << " Nphotons " << GetPhotonsNumber()
 434            << " Hough    " << nPhotonHough
 435            << " norm     " << nPhotonHoughNorm << endl;
 436
 437       cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl;
 438       cout << " ------------------------------------- " << endl;
 439     }
 440
 441   Int_t nPhotons = GetPhotonsNumber();
 442
 443   Float_t xmean = 0.;
 444   Float_t x2mean = 0.;
 445   Int_t nev = 0;
 446
 447   for (Int_t j=0; j < nPhotons;j++)
 448     {
 449       SetPhotonIndex(j);
 450
 451       Float_t eta = GetPhotonEta();
 452
 453       if(eta != -999.)
 454         {
 455           if(GetPhotonFlag() == 2)
 456             {
 457
 458
 459               xmean += eta;
 460               x2mean += eta*eta;
 461               nev++;
 462             }
 463         }
 464     }
 465
 466   if(nev > 0)
 467     {
 468       xmean /=(Float_t)nev;
 469       x2mean /=(Float_t)nev;
 470     } else {
 471       xmean = 0.;
 472       x2mean = 0.;
 473     }
 474
 475   Float_t vRMS = sqrt(x2mean - xmean*xmean);
 476
 477   SetHoughRMS(vRMS);
 478
 479   if(fDebug) cout << " RMS " << vRMS << endl;
 480
 481 }
 482
 483 void AliRICHRecon::FindEmissionPoint()
 484 {
 485   //estimate the emission point in radiator
 486
 487 // Find emission point
 488
 489   Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm)
 490   // 7.83 = -1/ln(T0) where
 491   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
 492   Float_t photonLenght, photonLenghtMin, photonLenghtMax;
 493
 494   photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
 495   photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght);
 496   photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
 497   Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
 498
 499   SetEmissionPoint(emissionPoint);
 500 }
 501
 502
 503 Int_t AliRICHRecon::PhotonInBand()
 504 {
 505   //search band fro photon candidates
 506
 507   //  Float_t massOfParticle;
 508   Float_t beta;
 509   Float_t nfreon;
 510
 511   Float_t thetacer;
 512
 513   Float_t xtoentr = GetEntranceX();
 514   Float_t ytoentr = GetEntranceY();
 515
 516   Float_t innerRadius;
 517   Float_t outerRadius;
 518
 519   Float_t phpad = GetPhiPoint();
 520
 521   //  Float_t pmod = GetTrackMomentum();
 522   //  Float_t trackTheta = GetTrackTheta();
 523   //  Float_t trackPhi = GetTrackPhi();
 524
 525   // inner radius //
 526   SetPhotonEnergy(5.6);
 527   SetEmissionPoint(fRadiatorWidth -0.0001);
 528   SetMassHypotesis(0.93828);
 529
 530   SetBetaOfParticle();
 531   SetFreonRefractiveIndex();
 532
 533   beta   = GetBetaOfParticle();
 534   nfreon = GetFreonRefractiveIndex();
 535
 536   thetacer = Cerenkovangle(nfreon,beta);
 537
 538   thetacer = 0.;
 539
 540   if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl;
 541
 542   FindThetaAtQuartz(thetacer);
 543
 544   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 545     {
 546       innerRadius = -999.;
 547       SetXInnerRing(-999.);
 548       SetYInnerRing(-999.);
 549       SetRadiusInnerRing(-999.);
 550     }
 551   else
 552     {
 553       SetThetaPhotonInDRS(GetThetaAtQuartz());
 554       SetPhiPhotonInDRS(phpad);
 555
 556       innerRadius = FromEmissionToCathode();
 557        if(innerRadius == 999.) innerRadius = -999.;
 558
 559       SetXInnerRing(GetXPointOnCathode());
 560       SetYInnerRing(GetYPointOnCathode());
 561       SetRadiusInnerRing(innerRadius);
 562     }
 563
 564   // outer radius //
 565   SetPhotonEnergy(7.7);
 566   SetEmissionPoint(0.);
 567 //  SetMassHypotesis(0.139567);
 568   SetMassHypotesis(0.);
 569
 570   SetBetaOfParticle();
 571   SetFreonRefractiveIndex();
 572
 573   beta   = GetBetaOfParticle();
 574   nfreon = GetFreonRefractiveIndex();
 575
 576   thetacer = Cerenkovangle(nfreon,beta);
 577
 578   //  thetacer = 0.75;
 579
 580   if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl;
 581
 582   FindThetaAtQuartz(thetacer);
 583
 584   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 585     {
 586       outerRadius = 999.;
 587       SetXOuterRing(999.);
 588       SetYOuterRing(999.);
 589       SetRadiusOuterRing(999.);
 590     }
 591   else
 592     {
 593       SetThetaPhotonInDRS(GetThetaAtQuartz());
 594       SetPhiPhotonInDRS(phpad);
 595
 596       outerRadius = FromEmissionToCathode();
 597 //      cout << " outerRadius " << outerRadius << endl;
 598       SetXOuterRing(GetXPointOnCathode());
 599       SetYOuterRing(GetYPointOnCathode());
 600       SetRadiusOuterRing(outerRadius);
 601     }
 602
 603   Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
 604
 605   if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius);
 606
 607   if(padradius>=innerRadius && padradius<=outerRadius) return 1;
 608   return 0;
 609 }
 610
 611 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
 612 {
 613   //find the theta at the quartz plate
 614
 615   if(thetaCerenkov == 999.)
 616     {
 617       SetThetaAtQuartz(999.);
 618       return;
 619     }
 620
 621   Float_t thetaAtQuartz = 999.;
 622
 623   Float_t trackTheta = GetTrackTheta();
 624
 625   if(trackTheta == 0) {
 626
 627     if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl;
 628
 629     thetaAtQuartz = thetaCerenkov;
 630     SetThetaAtQuartz(thetaAtQuartz);
 631     return;
 632   }
 633
 634   Float_t trackPhi   = GetTrackPhi();
 635   Float_t phiPoint = GetPhiPoint();
 636
 637   Double_t den = TMath::Sin((Double_t)trackTheta)
 638     *TMath::Cos((Double_t)trackPhi)
 639     *TMath::Cos((Double_t)phiPoint) +
 640     TMath::Sin((Double_t)trackTheta)
 641     *TMath::Sin((Double_t)trackPhi)
 642     *TMath::Sin((Double_t)phiPoint);
 643   Double_t b = TMath::Cos((Double_t)trackTheta)/den;
 644   Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den;
 645
 646   Double_t underSqrt = 1 + b*b - c*c;
 647
 648   if(fDebug)
 649     {
 650       cout << " trackTheta    " << trackTheta    << endl;
 651       cout << " TrackPhi      " << trackPhi      << endl;
 652       cout << " PhiPoint      " << phiPoint      << endl;
 653       cout << " ThetaCerenkov " << thetaCerenkov << endl;
 654       cout << " den b c " << den << " b " << b << " c " << c << endl;
 655     }
 656
 657   if(underSqrt < 0) {
 658     if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl;
 659     SetThetaAtQuartz(999.);
 660     return;
 661   }
 662
 663   Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
 664   Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
 665
 666   Double_t thetaSol1 = 2*TMath::ATan(sol1);
 667   Double_t thetaSol2 = 2*TMath::ATan(sol2);
 668
 669   if(fDebug) cout << " Theta sol 1 " << thetaSol1
 670                   << " Theta sol 2 " << thetaSol2 << endl;
 671
 672   if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
 673   if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
 674
 675   SetThetaAtQuartz(thetaAtQuartz);
 676 }
 677
 678 void AliRICHRecon::FindThetaPhotonCerenkov()
 679 {
 680   //find theta cerenkov of ring
 681
 682   Float_t thetaCerMin = 0.;
 683   Float_t thetaCerMax = 0.75;
 684   Float_t thetaCerMean;
 685
 686   Float_t radiusMin, radiusMax, radiusMean;
 687   Int_t nIteration = 0;
 688
 689   const Float_t kTollerance = 0.05;
 690
 691   //  Float_t pmod = GetTrackMomentum();
 692   //  Float_t trackTheta = GetTrackTheta();
 693   //  Float_t trackPhi = GetTrackPhi();
 694
 695   Float_t phiPoint = GetPhiPoint();
 696
 697   SetPhotonEnergy(6.85);
 698   SetEmissionPoint(fRadiatorWidth/2);
 699
 700   Float_t xPoint = GetEntranceX();
 701   Float_t yPoint = GetEntranceY();
 702   Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint);
 703
 704   if(fDebug) cout << " DistPoint " << distPoint << endl;
 705
 706   // Star minimization...
 707
 708   // First value...
 709
 710   FindThetaAtQuartz(thetaCerMin);
 711
 712   if(GetThetaAtQuartz() == 999.)
 713     {
 714       radiusMin = -999.;
 715     }
 716   else
 717     {
 718       SetThetaPhotonInDRS(GetThetaAtQuartz());
 719       SetPhiPhotonInDRS(phiPoint);
 720
 721       radiusMin = FromEmissionToCathode();
 722     }
 723
 724   // Second value...
 725
 726   FindThetaAtQuartz(thetaCerMax);
 727   if(GetThetaAtQuartz() == 999.)
 728     {
 729       radiusMax = 999.;
 730     }
 731   else
 732     {
 733       SetThetaPhotonInDRS(GetThetaAtQuartz());
 734       SetPhiPhotonInDRS(phiPoint);
 735
 736       radiusMax = FromEmissionToCathode();
 737     }
 738   // Mean value...
 739
 740   thetaCerMean = (thetaCerMax + thetaCerMin)/2;
 741
 742   FindThetaAtQuartz(thetaCerMean);
 743   if(GetThetaAtQuartz() == 999.)
 744     {
 745       radiusMean = 999.;
 746     }
 747   else
 748     {
 749       SetThetaPhotonInDRS(GetThetaAtQuartz());
 750       SetPhiPhotonInDRS(phiPoint);
 751
 752       radiusMean = FromEmissionToCathode();
 753     }
 754
 755   if(fDebug) cout << " r1 " << radiusMin << " rmean "
 756                   << radiusMean << " r2 " << radiusMax << endl;
 757
 758   while (TMath::Abs(radiusMean-distPoint) > kTollerance)
 759     {
 760
 761       if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
 762       if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
 763
 764         thetaCerMin = thetaCerMean;
 765
 766         FindThetaAtQuartz(thetaCerMin);
 767         SetThetaPhotonInDRS(GetThetaAtQuartz());
 768         SetPhiPhotonInDRS(phiPoint);
 769
 770         radiusMin =FromEmissionToCathode();
 771       }
 772
 773       thetaCerMean = (thetaCerMax + thetaCerMin)/2;
 774
 775       FindThetaAtQuartz(thetaCerMean);
 776       SetThetaPhotonInDRS(GetThetaAtQuartz());
 777       SetPhiPhotonInDRS(phiPoint);
 778
 779       radiusMean = FromEmissionToCathode();
 780
 781       nIteration++;
 782       if(nIteration>=50) {
 783         if(fDebug) printf(" max iterations in FindPhotonCerenkov\n");
 784         SetThetaPhotonCerenkov(999.);
 785         return;
 786       }
 787     }
 788
 789   SetThetaPhotonCerenkov(thetaCerMean);
 790
 791 }
 792
 793 void AliRICHRecon::FindAreaAndPortionOfRing()
 794 {
 795   //find fraction of the ring accepted by the RICH
 796
 797   Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
 798
 799   //  Float_t xtoentr = GetEntranceX();
 800   //  Float_t ytoentr = GetEntranceY();
 801   Float_t shiftX = GetShiftX();
 802   Float_t shiftY = GetShiftY();
 803
 804   Float_t xemiss = GetXCoordOfEmission();
 805   Float_t yemiss = GetYCoordOfEmission();
 806
 807   Float_t x0 = xemiss + shiftX;
 808   Float_t y0 = yemiss + shiftY;
 809
 810   //  Float_t pmod = GetTrackMomentum();
 811   //  Float_t trackTheta = GetTrackTheta();
 812   //  Float_t trackPhi = GetTrackPhi();
 813
 814   SetPhotonEnergy(6.85);
 815   SetFreonRefractiveIndex();
 816
 817   SetEmissionPoint(fRadiatorWidth/2.);
 818
 819   Float_t theta = GetThetaOfRing();
 820
 821   Int_t nPoints = 0;
 822   Int_t nPsiAccepted = 0;
 823   Int_t nPsiTotal = 0;
 824
 825   for(Int_t i=0;i<NPointsOfRing-1;i++)
 826     {
 827
 828       Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
 829
 830       SetThetaPhotonInTRS(theta);
 831       SetPhiPhotonInTRS(psi);
 832       FindPhotonAnglesInDRS();
 833
 834       Float_t radius = FromEmissionToCathode();
 835       if (radius == 999.) continue;
 836
 837       nPsiTotal++;
 838
 839       Float_t xPointRing = GetXPointOnCathode() + shiftX;
 840       Float_t yPointRing = GetYPointOnCathode() + shiftY;
 841
 842       SetDetectorWhereX(xPointRing);
 843       SetDetectorWhereY(yPointRing);
 844
 845       Int_t zone = CheckDetectorAcceptance();
 846
 847
 848       if (zone != 0)
 849         {
 850           FindIntersectionWithDetector();
 851           xPoint[nPoints] = GetIntersectionX();
 852           yPoint[nPoints] = GetIntersectionY();
 853         }
 854       else
 855         {
 856           xPoint[nPoints] = xPointRing;
 857           yPoint[nPoints] = yPointRing;
 858           nPsiAccepted++;
 859         }
 860
 861       nPoints++;
 862
 863     }
 864
 865   xPoint[nPoints] = xPoint[0];
 866   yPoint[nPoints] = yPoint[0];
 867
 868   // find area...
 869
 870   Float_t area = 0;
 871
 872   for (Int_t i = 0; i < nPoints; i++)
 873     {
 874       area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
 875     }
 876
 877   area *= 0.5;
 878
 879   Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
 880
 881
 882   SetAreaOfRing(area);
 883   SetPortionOfRing(portionOfRing);
 884 }
 885
 886 void AliRICHRecon::FindIntersectionWithDetector()
 887 {
 888   // find ring intersection with CsI edges
 889
 890   Float_t xIntersect, yIntersect;
 891   Float_t x1, x2, y1, y2;
 892
 893   Float_t shiftX = GetShiftX();
 894   Float_t shiftY = GetShiftY();
 895
 896   Float_t xPoint = GetXPointOnCathode() + shiftX;
 897   Float_t yPoint = GetYPointOnCathode() + shiftY;
 898
 899   Float_t xemiss = GetXCoordOfEmission();
 900   Float_t yemiss = GetYCoordOfEmission();
 901
 902   Float_t phi = GetPhiPhotonInDRS();
 903   Float_t m = tan(phi);
 904
 905   Float_t x0 = xemiss + shiftX;
 906   Float_t y0 = yemiss + shiftY;
 907
 908   if(xPoint > x0)
 909     {
 910       x1 = x0;
 911       x2 = xPoint;
 912     }
 913   else
 914     {
 915       x2 = x0;
 916       x1 = xPoint;
 917     }
 918   if(yPoint > y0)
 919     {
 920       y1 = y0;
 921       y2 = yPoint;
 922     }
 923   else
 924     {
 925       y2 = y0;
 926       y1 = yPoint;
 927     }
 928   //
 929   xIntersect = fXmax;
 930   yIntersect = m*(xIntersect - x0) + y0;
 931   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
 932     {
 933       SetIntersectionX(xIntersect);
 934       SetIntersectionY(yIntersect);
 935       return;
 936     }
 937   //
 938   xIntersect = fXmin;
 939   yIntersect = m*(xIntersect - x0) + y0;
 940   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
 941     {
 942       SetIntersectionX(xIntersect);
 943       SetIntersectionY(yIntersect);
 944       return;
 945     }
 946   //
 947   yIntersect = fYmax;
 948   xIntersect = (yIntersect - y0)/m + x0;
 949   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
 950     {
 951       SetIntersectionX(xIntersect);
 952       SetIntersectionY(yIntersect);
 953       return;
 954     }
 955   //
 956   yIntersect = fYmin;
 957   xIntersect = (yIntersect - y0)/m + x0;
 958   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
 959     {
 960       SetIntersectionX(xIntersect);
 961       SetIntersectionY(yIntersect);
 962       return;
 963     }
 964
 965   cout << " sono fuori!!!!!!" << endl;
 966
 967 }
 968 //__________________________________________________________________________________________________
 969 Int_t AliRICHRecon::CheckDetectorAcceptance() const
 970 {
 971   // check for the acceptance
 972
 973   // crosses X -2.6 2.6 cm
 974   // crosses Y -1 1 cm
 975
 976   Float_t xcoord = GetDetectorWhereX();
 977   Float_t ycoord = GetDetectorWhereY();
 978
 979   if(xcoord > fXmax)
 980     {
 981       if(ycoord > fYmax) return 2;
 982       if(ycoord > fYmin && ycoord < fYmax) return 3;
 983       if(ycoord < fYmin) return 4;
 984     }
 985   if(xcoord < fXmin)
 986     {
 987       if(ycoord > fYmax) return 8;
 988       if(ycoord > fYmin && ycoord < fYmax) return 7;
 989       if(ycoord < fYmin) return 6;
 990     }
 991   if(xcoord > fXmin && xcoord < fXmax)
 992     {
 993       if(ycoord > fYmax) return 1;
 994       if(ycoord > fYmin && ycoord < fYmax) return 0;
 995       if(ycoord < fYmin) return 5;
 996     }
 997   return 999;
 998 }
 999 //__________________________________________________________________________________________________
1000 Float_t AliRICHRecon::PhotonPositionOnCathode()
1001 {
1002   // find the photon position on the CsI
1003   //  Float_t massOfParticle;
1004   Float_t beta;
1005   Float_t nfreon;
1006
1007
1008   SetPhotonEnergy(6.85);
1009   SetEmissionPoint(fRadiatorWidth/2.);
1010   SetMassHypotesis(0.139567);
1011
1012   SetBetaOfParticle();
1013   SetFreonRefractiveIndex();
1014
1015   beta   = GetBetaOfParticle();
1016   nfreon = GetFreonRefractiveIndex();
1017
1018
1019   Float_t radius = FromEmissionToCathode();
1020   if (radius == 999.) return 999.;
1021
1022   return 0;
1023 }
1024
1025 void AliRICHRecon::FindPhotonAnglesInDRS()
1026 {
1027   // Setup the rotation matrix of the track...
1028
1029   TRotation mtheta;
1030   TRotation mphi;
1031   TRotation minv;
1032   TRotation mrot;
1033
1034   Float_t trackTheta = GetTrackTheta();
1035   Float_t trackPhi = GetTrackPhi();
1036
1037   mtheta.RotateY(trackTheta);
1038   mphi.RotateZ(trackPhi);
1039
1040   mrot = mphi * mtheta;
1041   //  minv = mrot.Inverse();
1042
1043   TVector3 photonInRadiator(1,1,1);
1044
1045   Float_t thetaCerenkov = GetThetaPhotonInTRS();
1046   Float_t phiCerenkov   = GetPhiPhotonInTRS();
1047
1048   photonInRadiator.SetTheta(thetaCerenkov);
1049   photonInRadiator.SetPhi(phiCerenkov);
1050   photonInRadiator = mrot * photonInRadiator;
1051   Float_t theta = photonInRadiator.Theta();
1052   Float_t phi = photonInRadiator.Phi();
1053   SetThetaPhotonInDRS(theta);
1054   SetPhiPhotonInDRS(phi);
1055
1056 }
1057
1058 Float_t AliRICHRecon::FromEmissionToCathode()
1059 {
1060   // trace from emission point to cathode
1061
1062   Float_t nfreon, nquartz, ngas;
1063
1064   SetFreonRefractiveIndex();
1065   SetQuartzRefractiveIndex();
1066   SetGasRefractiveIndex();
1067
1068   nfreon  = GetFreonRefractiveIndex();
1069   nquartz = GetQuartzRefractiveIndex();
1070   ngas    = GetGasRefractiveIndex();
1071
1072   Float_t trackTheta = GetTrackTheta();
1073   Float_t trackPhi = GetTrackPhi();
1074   Float_t lengthOfEmissionPoint = GetEmissionPoint();
1075
1076   Float_t theta = GetThetaPhotonInDRS();
1077   Float_t phi   = GetPhiPhotonInDRS();
1078
1079 //   cout << " Theta " << Theta << " Phi " << Phi << endl;
1080
1081   Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
1082   Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
1083
1084   SetXCoordOfEmission(xemiss);
1085   SetYCoordOfEmission(yemiss);
1086
1087   Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
1088
1089   if(thetaquar == 999.)
1090     {
1091       SetXPointOnCathode(999.);
1092       SetYPointOnCathode(999.);
1093       return thetaquar;
1094     }
1095
1096   Float_t thetagap  = SnellAngle( nquartz, ngas, thetaquar);
1097
1098   if(thetagap == 999.)
1099     {
1100       SetXPointOnCathode(999.);
1101       SetYPointOnCathode(999.);
1102       return thetagap;
1103     }
1104
1105   Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta);
1106   Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar);
1107   Float_t xg = fGapWidth*cos(phi)*tan(thetagap);
1108   Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta);
1109   Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar);
1110   Float_t yg = fGapWidth*sin(phi)*tan(thetagap);
1111
1112
1113   Float_t xtot = xemiss + xw + xq + xg;
1114   Float_t ytot = yemiss + yw + yq + yg;
1115
1116   SetXPointOnCathode(xtot);
1117   SetYPointOnCathode(ytot);
1118
1119
1120   Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
1121                                     +TMath::Power(fPhotonLimitY,2));
1122
1123   return distanceFromEntrance;
1124
1125 }
1126
1127
1128 void AliRICHRecon::FindPhiPoint()
1129 {
1130   //find phi of generated point
1131
1132   Float_t xtoentr = GetEntranceX();
1133   Float_t ytoentr = GetEntranceY();
1134
1135   Float_t trackTheta = GetTrackTheta();
1136   Float_t trackPhi = GetTrackPhi();
1137
1138   Float_t emissionPoint = GetEmissionPoint();
1139
1140   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
1141   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
1142   Float_t phipad = atan2(argY,argX);
1143
1144   SetPhiPoint(phipad);
1145
1146 }
1147
1148 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
1149 {
1150   // cerenkov angle from n and beta
1151
1152 // Compute the cerenkov angle
1153
1154   Float_t thetacer;
1155
1156   if((n*beta)<1.) {
1157     thetacer = 999.;
1158     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
1159     return thetacer;
1160   }
1161
1162   thetacer = acos (1./(n*beta));
1163   return thetacer;
1164 }
1165
1166 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
1167 {
1168   // Snell law
1169
1170 // Compute the Snell angle
1171
1172   Float_t sinrefractangle;
1173   Float_t refractangle;
1174
1175   sinrefractangle = (n1/n2)*sin(theta1);
1176
1177   if(sinrefractangle>1.) {
1178     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
1179     refractangle = 999.;
1180     return refractangle;
1181   }
1182
1183   refractangle = asin(sinrefractangle);
1184   return refractangle;
1185 }
1186
1187
1188 void AliRICHRecon::HoughResponse()
1189 {
1190   //Hough response
1191
1192 // Implement Hough response pat. rec. method
1193
1194   Float_t *hCSspace;
1195
1196   int           bin=0;
1197   int           bin1=0;
1198   int           bin2=0;
1199   int           i, j, k, nCorrBand;
1200   float         hcs[750],hcsw[750];
1201   float         angle, weight;
1202   float         lowerlimit,upperlimit;
1203
1204   float         etaPeak[100];
1205
1206   int           nBin;
1207
1208   float etaPeakPos  = -1;
1209
1210   Int_t   etaPeakCount = -1;
1211
1212   Float_t thetaCerenkov = 0.;
1213
1214   nBin = (int)(0.5+fThetaMax/(fDTheta));
1215   nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta));
1216
1217   memset ((void *)hcs, 0, fThetaBin*sizeof(float));
1218   memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
1219
1220   Int_t nPhotons = GetPhotonsNumber();
1221
1222   Int_t weightFlag = 0;
1223
1224   for (k=0; k< nPhotons; k++) {
1225
1226     SetPhotonIndex(k);
1227
1228     angle = GetPhotonEta();
1229
1230     if(angle == -999.) continue;
1231
1232     if (angle>=fThetaMin && angle<= fThetaMax)
1233
1234       {
1235
1236         bin = (int)(0.5+angle/(fDTheta));
1237
1238         bin1= bin-nCorrBand;
1239         bin2= bin+nCorrBand;
1240
1241         // calculate weights
1242
1243         if(fIsWEIGHT)
1244           {
1245             lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
1246             SetThetaOfRing(lowerlimit);
1247             FindAreaAndPortionOfRing();
1248             Float_t area1 = GetAreaOfRing();
1249
1250             upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
1251             SetThetaOfRing(upperlimit);
1252             FindAreaAndPortionOfRing();
1253             Float_t area2 = GetAreaOfRing();
1254
1255             //      cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
1256             Float_t diffarea = area2 - area1;
1257
1258             if(diffarea>0)
1259               {
1260                 weight = 1./(area2-area1);
1261               }
1262             else
1263               {
1264                 weightFlag = 1;
1265                 weight = 1.;
1266               }
1267
1268             //      cout <<" low "<< lowerlimit << " up " << upperlimit <<
1269             //        " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
1270
1271           }
1272         else
1273           {
1274             weight = 1.;
1275           }
1276
1277         SetPhotonWeight(weight);
1278
1279         //      cout << "weight..." << weight << endl;
1280
1281
1282         if (bin1<0)    bin1=0;
1283         if (bin2>nBin) bin2=nBin;
1284
1285         for (j=bin1; j<bin2; j++)
1286           {
1287             hcs[j] += 1;
1288             hcsw[j] += weight;
1289           }
1290       }
1291   }
1292
1293
1294   if(weightFlag == 0)
1295     {
1296       hCSspace = hcsw;
1297     }
1298   else
1299     {
1300       hCSspace = hcs;
1301       //      cout << " probems with weight...normal procedure adopted " << endl;
1302     }
1303
1304   HoughFiltering(hCSspace);
1305
1306   for (bin=0; bin <nBin; bin++) {
1307     angle = (bin+0.5) * (fDTheta);
1308     if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
1309       etaPeakCount = 0;
1310       etaPeakPos = hCSspace[bin];
1311       etaPeak[0]=angle;
1312     }
1313     else {
1314       if (hCSspace[bin] == etaPeakPos) {
1315         etaPeak[++etaPeakCount] = angle;
1316       }
1317     }
1318   }
1319
1320   for (i=0; i<etaPeakCount+1; i++) {
1321     thetaCerenkov += etaPeak[i];
1322   }
1323   if (etaPeakCount>=0) {
1324     thetaCerenkov /= etaPeakCount+1;
1325     fThetaPeakPos = etaPeakPos;
1326   }
1327
1328   SetThetaCerenkov(thetaCerenkov);
1329 }
1330
1331
1332 void AliRICHRecon::HoughFiltering(float hcs[])
1333 {
1334   // filter for Hough
1335
1336 // hough filtering
1337
1338    float hcsFilt[750];
1339    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
1340    int nx, i, nxDx;
1341    int sizeHCS;
1342    int nBin;
1343
1344    nBin =  (int)(1+fThetaMax/fDTheta);
1345    sizeHCS = fThetaBin*sizeof(float);
1346
1347    memset ((void *)hcsFilt, 0, sizeHCS);
1348
1349    for (nx = 0; nx < nBin; nx++) {
1350       for (i = 0; i < 5; i++)   {
1351         nxDx = nx + (i-2);
1352         if (nxDx> -1 && nxDx<nBin)
1353              hcsFilt[nx] +=  hcs[nxDx] * k[i];
1354       }
1355    }
1356
1357    for (nx = 0; nx < nBin; nx++) {
1358      hcs[nx] = hcsFilt[nx];
1359    }
1360 }
1361
1362 void AliRICHRecon::FindWeightThetaCerenkov()
1363 {
1364   // manage with weight for photons
1365
1366   Float_t wei = 0.;
1367   Float_t weightThetaCerenkov = 0.;
1368
1369   Int_t nPhotons = GetPhotonsNumber();
1370   for(Int_t i=0;i<nPhotons;i++)
1371     {
1372       SetPhotonIndex(i);
1373
1374       if(GetPhotonFlag() == 2)
1375         {
1376           Float_t photonEta = GetPhotonEta();
1377           Float_t photonWeight = GetPhotonWeight();
1378           weightThetaCerenkov += photonEta*photonWeight;
1379           wei += photonWeight;
1380         }
1381     }
1382
1383   if(wei != 0.)
1384     {
1385       weightThetaCerenkov /= wei;
1386     }
1387   else
1388     {
1389       weightThetaCerenkov = 0.;
1390     }
1391
1392   SetThetaCerenkov(weightThetaCerenkov);
1393
1394   cout << " thetac weighted -> " << weightThetaCerenkov << endl;
1395 }
1396
1397
1398 void AliRICHRecon::FlagPhotons()
1399 {
1400   // flag photons
1401
1402   Int_t nPhotonHough = 0;
1403
1404   Float_t thetaCerenkov = GetThetaCerenkov();
1405   if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl;
1406
1407   Float_t thetaDist= thetaCerenkov - fThetaMin;
1408   Int_t steps = (Int_t)(thetaDist / fDTheta);
1409
1410   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
1411   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
1412   Float_t tavg = 0.5*(tmin+tmax);
1413
1414   tmin = tavg - 0.5*fWindowWidth;
1415   tmax = tavg + 0.5*fWindowWidth;
1416
1417   if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl;
1418   if(fDebug) cout << " thetac " << thetaCerenkov << endl;
1419
1420   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
1421
1422   Int_t nPhotons = GetPhotonsNumber();
1423
1424   //  for(Int_t i=0;i<candidatePhotonsNumber;i++)
1425
1426   for(Int_t i=0;i<nPhotons;i++)
1427     {
1428       SetPhotonIndex(i);
1429
1430       Float_t photonEta = GetPhotonEta();
1431
1432       if(photonEta == -999.) continue;
1433
1434       if(photonEta >= tmin && photonEta <= tmax)
1435         {
1436           SetPhotonFlag(2);
1437           nPhotonHough++;
1438         }
1439     }
1440   SetHoughPhotons(nPhotonHough);
1441 }
1442
1443 void AliRICHRecon::DrawEvent(Int_t flag) const
1444 {
1445   // draw event with rings
1446
1447   flag=1; // dummy to be removed...
1448 }
1449
1450 Float_t  AliRICHRecon::FindMassOfParticle()
1451 {
1452   // find mass of the particle from theta cerenkov
1453
1454   Float_t pmod = GetTrackMomentum();
1455
1456   SetPhotonEnergy(6.85);
1457   SetFreonRefractiveIndex();
1458
1459   Float_t thetaCerenkov = GetThetaCerenkov();
1460   FindBetaFromTheta(thetaCerenkov);
1461
1462   Double_t beta = (Double_t)(GetBetaOfParticle());
1463   Double_t den = 1. - beta*beta;
1464   if(den<=0.) return 999.;
1465
1466   Double_t gamma = 1./TMath::Sqrt(den);
1467
1468   Float_t mass = pmod/(beta*(Float_t)gamma);
1469
1470   return mass;
1471 }
1472
1473
1474 void AliRICHRecon::FillHistograms()
1475 {
1476   // fill histograms..
1477
1478   Float_t fittedTrackTheta, fittedTrackPhi;
1479
1480   Float_t thetaCerenkov    = GetThetaCerenkov();
1481   if(thetaCerenkov == 999.) return;
1482
1483   Float_t vertZ = GetEventVertexZ();
1484
1485   Float_t trackTheta = GetTrackTheta();
1486   Float_t trackPhi   = GetTrackPhi();
1487   Float_t pmod       = GetTrackMomentum();
1488   Float_t pt         = GetTrackPt();
1489   Float_t trackEta   = GetTrackEta();
1490   Int_t q            = GetTrackCharge();
1491   Float_t tPCLastZ   = GetTrackTPCLastZ();
1492   Float_t minDist    = GetMinDist();
1493
1494   fittedTrackTheta = GetFittedTrackTheta();
1495   fittedTrackPhi   = GetFittedTrackPhi();
1496   Int_t fittednPhotonHough = GetFittedHoughPhotons();
1497
1498   if(fDebug)
1499     {
1500       cout << " p " << pmod  << " ThetaC " << thetaCerenkov
1501            << " rings " << fNrings << endl;
1502     }
1503
1504   Int_t nPhotonHough     = GetHoughPhotons();
1505 //  Float_t nPhotonHoughNorm = GetHoughPhotonsNorm();
1506   Float_t inRing = GetPortionOfRing();
1507
1508   Float_t massOfParticle = FindMassOfParticle();
1509
1510   Float_t houghArea = GetHoughArea();
1511   Float_t multiplicity = GetEventMultiplicity();
1512
1513
1514   Float_t var[20];
1515
1516   var[0] = 0;
1517   var[1] = 0;
1518   var[2] = vertZ;
1519   var[3] = pmod;
1520   var[4] = pt;
1521   var[5] = trackEta;
1522   var[6] = trackTheta;
1523   var[7] = trackPhi;
1524   var[8] = fittedTrackTheta;
1525   var[9] = fittedTrackPhi;
1526   var[10] = q;
1527   var[11] = thetaCerenkov;
1528   var[12] = (Float_t)nPhotonHough;
1529   var[13] = (Float_t)fittednPhotonHough;
1530   var[14] = inRing;
1531   var[15] = massOfParticle;
1532   var[16] = houghArea;
1533   var[17] = multiplicity;
1534   var[18] = tPCLastZ;
1535   var[19] = minDist;
1536
1537   fNtuple->Fill(var);
1538
1539
1540   fittedTrackTheta = GetFittedTrackTheta();
1541   fittedTrackPhi = GetFittedTrackPhi();
1542
1543
1544
1545   if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) {
1546       SetPhotonEnergy(6.85);
1547       SetFreonRefractiveIndex();
1548   }
1549
1550   Int_t nPhotons = GetPhotonsNumber();
1551
1552   for (Int_t j=0; j < nPhotons;j++)
1553     SetPhotonIndex(j);
1554 }//FillHistograms()
1555 //__________________________________________________________________________________________________
1556 void AliRICHRecon::Minimization()
1557 {
1558   // minimization to find the best theta and phi of the track
1559
1560   Double_t arglist;
1561   Int_t ierflag = 0;
1562
1563   static Double_t vstart[2];
1564   static Double_t lower[2], upper[2];
1565   static Double_t step[2]={0.001,0.001};
1566
1567   Double_t trackThetaNew,trackPhiNew;
1568   TString chname;
1569   Double_t eps, b1, b2;
1570   Int_t ierflg;
1571
1572   gAliRICHminuit = new TMinuit(2);
1573   gAliRICHminuit->SetObjectFit((TObject *)this);
1574   gAliRICHminuit->SetFCN(fcnrecon);
1575   gAliRICHminuit->mninit(5,10,7);
1576
1577   vstart[0] = (Double_t)GetTrackTheta();
1578   vstart[1] = (Double_t)GetTrackPhi();
1579
1580   lower[0] = vstart[0] - 0.03;
1581   if(lower[0] < 0) lower[0] = 0.;
1582   upper[0] = vstart[0] + 0.03;
1583   lower[1] = vstart[1] - 0.03;
1584   upper[1] = vstart[1] + 0.03;
1585
1586
1587   gAliRICHminuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag);
1588   gAliRICHminuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag);
1589
1590   arglist = -1;
1591
1592   //  gAliRICHminuit->FixParameter(0);
1593
1594   gAliRICHminuit->SetPrintLevel(-1);
1595 //  gAliRICHminuit->mnexcm("SET PRI",&arglist, 1, ierflag);
1596   gAliRICHminuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
1597   gAliRICHminuit->mnexcm("SET NOW",&arglist, 1, ierflag);
1598   arglist = 1;
1599   gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg);
1600   arglist = -1;
1601
1602   //  gAliRICHminuit->mnscan();
1603
1604 //  gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
1605   gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
1606   gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
1607
1608   gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg);
1609   gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg);
1610
1611   //values after the fit...
1612   SetFittedTrackTheta((Float_t)trackThetaNew);
1613   SetFittedTrackPhi((Float_t)trackPhiNew);
1614
1615   delete gAliRICHminuit;
1616
1617 }
1618
1619 void AliRICHRecon::EstimationOfTheta()
1620 {
1621   // theta estimate
1622
1623   Int_t nPhotons = 0;
1624
1625   Float_t shiftX = GetShiftX();
1626   Float_t shiftY = GetShiftY();
1627
1628   Float_t *candidatePhotonX = GetCandidatePhotonX();
1629   Float_t *candidatePhotonY = GetCandidatePhotonY();
1630
1631   Int_t nPhotonsCandidates = GetCandidatePhotonsNumber();
1632
1633   //  cout << "MINIM: Nphotons " << nPhotonsCandidates << endl;
1634
1635   for (Int_t j=0; j < nPhotonsCandidates; j++)
1636     {
1637
1638       SetPhotonIndex(j);
1639
1640       if(!GetPhotonFlag()) continue;
1641
1642       Float_t xtoentr = candidatePhotonX[j] - shiftX;
1643       Float_t ytoentr = candidatePhotonY[j] - shiftY;
1644
1645       SetEntranceX(xtoentr);
1646       SetEntranceY(ytoentr);
1647
1648       FindPhiPoint();
1649
1650       FindThetaPhotonCerenkov();
1651
1652       Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
1653
1654       //      cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl;
1655
1656       SetPhotonEta(thetaPhotonCerenkov);
1657
1658       nPhotons++;
1659
1660     }
1661
1662   Float_t xmean = 0.;
1663   Float_t x2mean = 0.;
1664   Int_t nev = 0;
1665
1666   for (Int_t j=0; j < nPhotonsCandidates;j++)
1667     {
1668       SetPhotonIndex(j);
1669
1670       Float_t eta = GetPhotonEta();
1671
1672       if(eta != -999.)
1673         {
1674           if(GetPhotonFlag() == 2)
1675             {
1676               xmean += eta;
1677               x2mean += eta*eta;
1678               nev++;
1679             }
1680         }
1681     }
1682
1683   if(nev > 0)
1684     {
1685       xmean /=(Float_t)nev;
1686       x2mean /=(Float_t)nev;
1687     } else {
1688       xmean = 0.;
1689       x2mean = 0.;
1690     }
1691
1692   Float_t vRMS = sqrt(x2mean - xmean*xmean);
1693
1694   //  cout << " RMS " << vRMS;
1695
1696   SetEstimationOfTheta(xmean);
1697   SetEstimationOfThetaRMS(vRMS);
1698 }
1699
1700 void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t)
1701 {
1702   // function to be minimized
1703   AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit();
1704
1705   Float_t p0 = (Float_t)par[0];
1706   Float_t p1 = (Float_t)par[1];
1707
1708   gMyRecon->SetTrackTheta(p0);
1709   gMyRecon->SetTrackPhi(p1);
1710
1711   gMyRecon->EstimationOfTheta();
1712   Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS();
1713
1714   Int_t houghPhotons = gMyRecon->GetHoughPhotons();
1715
1716
1717   f = (Double_t)(1000*vRMS/(Float_t)houghPhotons);
1718
1719 //   if(fDebug) cout << "   f   " << f
1720 //                << " theta " << par[0] << " phi " << par[1]
1721 //                   << " HoughPhotons " << houghPhotons << endl;
1722 //
1723 //   if(fDebug&&iflag == 3)
1724 //     {
1725 //             cout << " --- end convergence...summary --- " << endl;
1726 //             cout << " theta " << par[0] << endl;
1727 //             cout << "  phi  " << par[1] << endl;
1728 //     }
1729 }
1730
1731 void AliRICHRecon::Waiting()
1732 {
1733   // wait, wait....
1734   if(!fIsDISPLAY) return;
1735   cout << " Press any key to continue...";
1736
1737 //  gSystem->ProcessEvents();
1738   getchar();
1739
1740   cout << endl;
1741
1742   return;
1743 }
1744