RICH/AliRICHRecon.cxx

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