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