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