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 #include "AliRICH.h"
  17 #include "AliRICHRecon.h"
  18 #include "AliRICHParam.h"
  19 #include <AliLoader.h>
  20 #include <AliStack.h>
  21 #include <Riostream.h>
  22 #include <TCanvas.h>
  23 #include <TParticle.h>
  24 #include <TStyle.h>
  25 #include <TF1.h>
  26 #include <TH2.h>
  27 #include <TMath.h>
  28 #include <TRandom.h>
  29 #include <TMinuit.h>
  30 #include <TNtuple.h>
  31 #include <TMath.h>
  32 #include <TGraph.h>
  33 #include <TLine.h>
  34 #include <TPolyLine.h>
  35 #include <TMarker.h>
  36 #include <TText.h>
  37 #include <TProfile.h>
  38 #include <TRotation.h>
  39 #include <TSystem.h>
  40 #include <TVector3.h>
  41 #include <TEventList.h>
  42
  43 #define NPointsOfRing 201
  44
  45 // Geometry of the RICH at Star...
  46
  47 static const Int_t nPadX      = AliRICHParam::NpadsY();
  48 static const Int_t nPadY      = AliRICHParam::NpadsX();
  49 static const Float_t PadSizeX = AliRICHParam::PadSizeY();
  50 static const Float_t PadSizeY = AliRICHParam::PadSizeX();
  51 static const Float_t spacer   = AliRICHParam::DeadZone();
  52 static const Float_t degree = 180/3.1415926535;
  53
  54 static const Float_t pi = TMath::Pi();
  55
  56 static const Float_t RadiatorWidth = AliRICHParam::FreonThickness();
  57 static const Float_t QuartzWidth   = AliRICHParam::QuartzThickness();
  58 static const Float_t GapWidth      = AliRICHParam::RadiatorToPads();
  59
  60 static const Float_t fDTheta = 0.001;          // Step for sliding window
  61 //static const Float_t fWindowWidth = 0.040;     // Hough width of sliding window
  62 static const Float_t fWindowWidth = 0.060;     // Hough width of sliding window
  63
  64 static const Int_t fThetaBin = 750;            // Hough bins
  65 static const Float_t fThetaMin = 0.0;          // Theta band min
  66 static const Float_t fThetaMax = 0.75;         // Theta band max
  67
  68 static const Float_t Xmin = -AliRICHParam::PcSizeY()/2.;
  69 static const Float_t Xmax =  AliRICHParam::PcSizeY()/2.;
  70 static const Float_t Ymin = -AliRICHParam::PcSizeX()/2.;
  71 static const Float_t Ymax =  AliRICHParam::PcSizeX()/2.;
  72
  73
  74 // Global variables...
  75
  76 Bool_t fDebug      = kFALSE;
  77 Bool_t kDISPLAY    = kFALSE;
  78 Bool_t kWEIGHT     = kFALSE;
  79 Bool_t kBACKGROUND = kFALSE;
  80 Bool_t kMINIMIZER  = kFALSE;
  81 //
  82
  83 Int_t TotEvents = 0;
  84
  85 static Float_t xGraph[3000],yGraph[3000];
  86
  87 static Int_t NRings = 0;
  88 static Int_t NevTOT = 0;
  89
  90 TMinuit *gMyMinuit ;
  91
  92 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  93
  94 // Float_t fEmissionPoint;
  95 // Float_t fTrackTheta;
  96 // Float_t fTrackPhi;
  97 // Float_t fXtoentr;
  98 // Float_t fYtoentr;
  99
 100 //
 101
 102 TFile *outputfile;
 103
 104 TH1F *h1_photons,*h1_photacc,*h1_hough;
 105 TH2F *h2_tvsppos, *h2_tvspneg,*h2_func;
 106
 107 TH2F *h2_disp;
 108
 109 TH2F *h2_test1, *h2_test2, *h2_test4, *h2_testmap;
 110 TH2F *h2_dist_p;
 111
 112 TH1F *h1_photons1, *h1_photons2;
 113 TH1F *h1_houghpos, *h1_houghneg;
 114 TH1F *h1_mass;
 115 TH2F *h2_mvsp;
 116
 117 TH1F *h1_hcs, *h1_hcsw;
 118
 119 TH1F *h1_nprotons;
 120
 121 TProfile *hp_1pos, *hp_1neg;
 122 TProfile *hp_1posnorm, *hp_1negnorm;
 123 TH2F *h2_1pos, *h2_1neg;
 124 TH2F *h2_1posnorm, *h2_1negnorm;
 125 TH2F *h2_mvst;
 126
 127 TH1F *h1_deltap, *h1_deltapop;
 128 TH1F *h1_diffTrackTheta, *h1_diffTrackPhi;
 129 TH1F *h1_photaccspread;
 130
 131 TH2F *h2_diffpos, *h2_diffneg;
 132 TH2F *h2_map, *h2_mapw;
 133
 134 TH1F *photris;
 135
 136 TNtuple *hn;
 137
 138 TCanvas *StarCanvas,*Display,*Displayhcs;
 139 TGraph *gra;
 140 TLine *line;
 141 TPolyLine *poll;
 142 TPolyMarker *polm;
 143 TMarker *Point, *TrackPoints, *Photon, *PhotonAcc;
 144 TText *text;
 145
 146 AliRICHRecon::AliRICHRecon(const char*, const char*)
 147 {
 148
 149   fRich = (AliRICH*)gAlice->GetDetector("RICH");
 150
 151 }
 152
 153 void AliRICHRecon::InitRecon()
 154 {
 155
 156    outputfile = new TFile("Anal.root","RECREATE","My Analysis histos");
 157    if(kDISPLAY) Display = new TCanvas("Display","RICH Display",0,0,1200,750);
 158
 159    h1_photons = new TH1F("h1_photons","photons",750,0.,0.75);
 160    h1_photacc = new TH1F("h1_photacc","photons",750,0.,0.75);
 161    h1_hough   = new TH1F("h1_hough","hough",750,0.,0.75);
 162    h1_houghpos= new TH1F("h1_houghpos","hough",750,0.,0.75);
 163    h1_houghneg= new TH1F("h1_houghneg","hough",750,0.,0.75);
 164
 165    h2_tvsppos = new TH2F("h2_tvsppos","thetac vs p",100,0.,5.,750,0.,0.75);
 166    h2_tvspneg = new TH2F("h2_tvspneg","thetac vs p",100,0.,5.,750,0.,0.75);
 167    h2_func    = new TH2F("h2_func"," func ",800,0.,0.8,100,-100.,100.);
 168    h2_mvsp    = new TH2F("h2_mvsp","mass vs p",100,0.,5.,200,0.,2.);
 169    h2_mvst    = new TH2F("h2_mvst","mass vs t",750,0.,0.75,200,0.,2.);
 170    h2_map     = new TH2F("h2_map","h2_map",160,0.,160.,96,0.,96.);
 171    h2_mapw    = new TH2F("h2_mapw","h2_mapw",160,0.,160.,96,0.,96.);
 172
 173    h2_dist_p = new TH2F("h2_dist_p","h2_dist_p",100,0.,5.,100,0.,5.);
 174    //
 175
 176    h2_disp = new TH2F("h2_disp","STAR-RICH Event Display",165,Xmin,Xmax,100,Ymin,Ymax);
 177
 178    //   h2_test1  = new TH2F("h2_test1","test1 map",165,-64.,64.,100,-42.,42.);
 179    h2_test2  = new TH2F("h2_test2","test2 map",165,-64.,64.,100,-42.,42.);
 180    //   h2_test4  = new TH2F("h2_test4","test4 map",165,-64.,64.,100,-42.,42.);
 181    h2_testmap= new TH2F("h2_testmap","test map",165,-64.,64.,100,-42.,42.);
 182
 183    //
 184    h1_photons1 = new TH1F("h1_photons1","photons",750,0.,0.75);
 185    h1_photons2 = new TH1F("h1_photons2","photons",750,0.,0.75);
 186    //
 187    h1_hcs  = new TH1F("h1_hcs","hcs",750,0.,750.);
 188    h1_hcsw = new TH1F("h1_hcsw","hcsw",750,0.,750.);
 189    //
 190    h1_nprotons = new TH1F("h1_nprotons","n prot",30,0.,30.);
 191    //
 192    hp_1pos = new TProfile("hp_1pos","Nphot vs thetac pos",250,0.,0.75);
 193    hp_1neg = new TProfile("hp_1neg","Nphot vs thetac neg",250,0.,0.75);
 194    hp_1posnorm = new TProfile("hp_1posnorm","Nphot vs thetac pos norm",250,0.,0.75);
 195    hp_1negnorm = new TProfile("hp_1negnorm","Nphot vs thetac neg norm",250,0.,0.75);
 196    //
 197    h2_1pos     = new TH2F("h2_1pos","Nphot vs p pos",100,0.,5.,30,0.,30.);
 198    h2_1neg     = new TH2F("h2_1neg","Nphot vs p neg",100,0.,5.,30,0.,30.);
 199    h2_1posnorm = new TH2F("h2_1posnorm","Nphot vs p pos norm",100,0.,5.,30,0.,30.);
 200    h2_1negnorm = new TH2F("h2_1negnorm","Nphot vs p neg norm",100,0.,5.,30,0.,30.);
 201
 202    h1_deltap = new TH1F("h1_deltap","delta_p",200,-0.5,0.5);
 203    h1_deltapop = new TH1F("h1_deltapop","deltapop",200,-1.,1.);
 204    h1_diffTrackTheta = new TH1F("h1_diffTrackTheta","delta theta",200,-0.25,0.25);
 205    h1_diffTrackPhi   = new TH1F("h1_diffTrackPhi","delta phi",200,-0.25,0.25);
 206
 207    h1_photaccspread = new TH1F("h1_photaccspread","photons spread",200,-0.1,0.1);
 208
 209    //
 210
 211    h1_mass = new TH1F("h1_mass","mass",200,0.,2.);
 212    photris = new TH1F("photris","photris",1000,0.,1.);
 213    h2_diffneg = new TH2F("h2_diffneg","diff neg",100,-2.5,2.5,100,-2.5,2.5);
 214    h2_diffpos = new TH2F("h2_diffpos","diff pos",100,-2.5,2.5,100,-2.5,2.5);
 215
 216    hn = new TNtuple("hn","ntuple",
 217 "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge:ThetaCerenkov:NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ");
 218 }
 219
 220 void AliRICHRecon::StartProcessEvent()
 221 {
 222
 223   Float_t TrackThetaStored    = 0;
 224   Float_t TrackPhiStored      = 0;
 225   Float_t ThetaCerenkovStored = 0;
 226   Int_t HoughPhotonsStored    = 0;
 227
 228   SetFreonScaleFactor(0.994);
 229
 230   InitRecon();
 231
 232   if(kDISPLAY)
 233     {
 234       DrawEvent(0);
 235 //      waiting();
 236     }
 237
 238     Rich()->GetLoader()->LoadHits();
 239     Rich()->GetLoader()->LoadRecPoints();
 240     Rich()->GetLoader()->LoadDigits();
 241     gAlice->GetRunLoader()->LoadHeader();
 242     gAlice->GetRunLoader()->LoadKinematics();
 243
 244     Rich()->GetLoader()->TreeR()->GetEntry(0);
 245
 246     Float_t clusX[7][500],clusY[7][500];
 247     Int_t clusQ[7][500];
 248     Int_t nClusters[7];
 249
 250     for (Int_t ich=0;ich<7;ich++) {
 251       nClusters[ich] = Rich()->ClustersOld(ich+1)->GetEntries();
 252       for(Int_t k=0;k<nClusters[ich];k++) {
 253         AliRICHRawCluster *pCluster = (AliRICHRawCluster *)Rich()->ClustersOld(ich+1)->At(k);
 254         clusX[ich][k] = pCluster->fX;
 255         clusY[ich][k] = pCluster->fY;
 256         clusQ[ich][k] = pCluster->fQ;
 257 //        pCluster->Print();
 258       }
 259     }
 260
 261     Int_t nPrimaries = (Int_t)Rich()->GetLoader()->TreeH()->GetEntries();
 262
 263     cout << " N. primaries " << nPrimaries << endl;
 264
 265     for(Int_t i=0;i<nPrimaries;i++){
 266
 267       Rich()->GetLoader()->TreeH()->GetEntry(i);
 268
 269       Rich()->Hits()->Print();
 270       Int_t iPrim = 0;
 271
 272       AliRICHhit* pHit=0;
 273
 274       for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) {
 275
 276         pHit = (AliRICHhit*)Rich()->Hits()->At(j);
 277         if(pHit->GetTrack() < nPrimaries) break;
 278         iPrim++;
 279       }
 280
 281       cout << " iPrim " << iPrim << endl;
 282 //      if(iPrim==0) return;
 283 //      if(iPrim>1) Fatal("StartProcessEvent"," problems with prim to hit!!! = %3i", iPrim);
 284
 285       TParticle *pParticle = gAlice->GetRunLoader()->Stack()->Particle(pHit->GetTrack());
 286       Float_t pmod     = pParticle->P();
 287       Float_t pt       = pParticle->Pt();
 288       Float_t TrackEta = pParticle->Eta();
 289       Int_t q          = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge());
 290
 291       pParticle->Print();
 292
 293       cout << " pmod " << pmod << " pt " << pt << " Eta " << TrackEta << " charge " << q << endl;
 294
 295       SetTrackMomentum(pmod);
 296       SetTrackPt(pt);
 297       SetTrackEta(TrackEta);
 298       SetTrackCharge(q);
 299
 300       TVector3 pGlob(pHit->MomFreoX(),pHit->MomFreoY(),pHit->MomFreoZ());
 301       TVector3 pLocal = Rich()->C(pHit->Chamber())->Global2Local(pGlob,1);
 302
 303       Float_t primGlobalX = pHit->X();
 304       Float_t primGlobalY = pHit->Y();
 305       Float_t primGlobalZ = pHit->Z();
 306       TVector3 primGlobal(primGlobalX,primGlobalY,primGlobalZ);
 307       TVector3 primLocal = Rich()->C(pHit->Chamber())->Global2Local(primGlobal);
 308
 309 //      Float_t pmodFreo = pLocal.Mag();
 310       Float_t TrackTheta = pLocal.Theta();
 311       Float_t TrackPhi = pLocal.Phi();
 312
 313       cout << " TrackTheta " << TrackTheta << " TrackPhi " << TrackPhi << endl;
 314
 315       SetTrackTheta(TrackTheta);
 316       SetTrackPhi(TrackPhi);
 317
 318       Int_t MaxInd = 0;
 319       Float_t MinDist =  999.;
 320
 321       cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl;
 322
 323       for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++)
 324         {
 325           Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j];
 326           Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j];
 327
 328           cout << " cluster x " << clusX[pHit->Chamber()-1][j] << " hit track x " << primLocal.X();
 329           cout << " cluster y " << clusY[pHit->Chamber()-1][j] << " hit track y " << primLocal.Y() << endl;
 330
 331           Float_t diff = sqrt(diffx*diffx + diffy*diffy);
 332
 333           if(diff < MinDist)
 334             {
 335               MinDist = diff;
 336               MaxInd = j;
 337             }
 338
 339         }
 340
 341       Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][MaxInd];
 342       Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][MaxInd];
 343
 344       cout << " diffx " << diffx << " diffy " << diffy << endl;
 345
 346       if(q>0)
 347       {
 348          h2_diffpos->Fill(diffx,diffy);
 349       } else {
 350          h2_diffneg->Fill(diffx,diffy);
 351       }
 352
 353       SetMipIndex(MaxInd);
 354       SetTrackIndex(i);
 355
 356       Float_t ShiftX = primLocal.X()/primLocal.Z()*(RadiatorWidth+QuartzWidth+GapWidth) + primLocal.X();
 357       Float_t ShiftY = primLocal.Y()/primLocal.Z()*(RadiatorWidth+QuartzWidth+GapWidth) + primLocal.Y();
 358
 359       SetShiftX(ShiftX);
 360       SetShiftY(ShiftY);
 361
 362       Float_t *pclusX = &clusX[pHit->Chamber()-1][0];
 363       Float_t *pclusY = &clusY[pHit->Chamber()-1][0];
 364
 365       SetCandidatePhotonX(pclusX);
 366       SetCandidatePhotonY(pclusY);
 367       SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]);
 368
 369       Int_t qch = clusQ[pHit->Chamber()-1][MaxInd];
 370
 371       if(MinDist < 3.0 && qch > 120 && MaxInd !=0)
 372         {
 373
 374           if(kBACKGROUND)
 375             {
 376
 377               Float_t Xrndm = Xmin + (Xmax-Xmin)*gRandom->Rndm(280964);
 378               Float_t Yrndm = Ymin + (Ymax-Ymin)*gRandom->Rndm(280964);
 379
 380               cout << " Xrndm " << Xrndm << " Yrndm " << Yrndm << endl;
 381
 382               SetShiftX(Xrndm);
 383               SetShiftY(Yrndm);
 384
 385             }
 386
 387           PatRec();
 388
 389           TrackThetaStored = GetTrackTheta();
 390           TrackPhiStored = GetTrackPhi();
 391           ThetaCerenkovStored = GetThetaCerenkov();
 392           HoughPhotonsStored = GetHoughPhotons();
 393
 394           Int_t DiffNPhotons = 999;
 395           Int_t Nsteps = 0;
 396           Float_t DiffTrackTheta = 999.;
 397           Float_t DiffTrackPhi   = 999.;
 398
 399           while( kMINIMIZER && GetHoughPhotons() > 2
 400                             && DiffNPhotons !=0
 401                             && DiffTrackTheta > 0.0001
 402                             && Nsteps < 10)
 403             {
 404
 405               Int_t   HoughPhotonsBefore  = GetHoughPhotons();
 406
 407               Float_t TrackThetaBefore = GetTrackTheta();
 408               Float_t TrackPhiBefore   = GetTrackPhi();
 409
 410               Minimization();
 411
 412               PatRec();
 413
 414               DiffNPhotons = abs(HoughPhotonsBefore - GetHoughPhotons());
 415
 416               Float_t TrackThetaAfter = GetTrackTheta();
 417               Float_t TrackPhiAfter   = GetTrackPhi();
 418
 419               DiffTrackTheta = TMath::Abs(TrackThetaAfter - TrackThetaBefore);
 420               DiffTrackPhi   = TMath::Abs(TrackPhiAfter - TrackPhiBefore);
 421
 422               if(fDebug)
 423               cout << " HoughPhotonsBefore " << HoughPhotonsBefore
 424                    << " GetHoughPhotons()  " << GetHoughPhotons();
 425
 426               Nsteps++;
 427             }
 428
 429           SetFittedThetaCerenkov(GetThetaCerenkov());
 430           SetFittedHoughPhotons(GetHoughPhotons());
 431
 432           SetTrackTheta(TrackThetaStored);
 433           SetTrackPhi(TrackPhiStored);
 434           SetThetaCerenkov(ThetaCerenkovStored);
 435           SetHoughPhotons(HoughPhotonsStored);
 436
 437           SetMinDist(MinDist);
 438
 439           FillHistograms();
 440
 441           if(kDISPLAY) DrawEvent(1);
 442
 443           waiting();
 444
 445         }
 446     }
 447   //
 448   if(kDISPLAY) Display->Print("display.ps");
 449 }
 450
 451
 452 void AliRICHRecon::EndProcessEvent()
 453 {
 454 // function called at the end of the event loop
 455
 456   printf("Processed events: %d Total events: %d \n",TotEvents,NevTOT);
 457
 458   outputfile->Write();
 459   outputfile->Close();
 460 }
 461
 462 void AliRICHRecon::PatRec()
 463 {
 464
 465   Float_t TrackTheta = GetTrackTheta();
 466   Float_t TrackPhi   = GetTrackPhi();
 467   Float_t pmod       = GetTrackMomentum();
 468   //  Int_t q            = GetTrackCharge();
 469
 470   //  Int_t TrackIndex = GetTrackIndex();
 471   Int_t MipIndex   = GetMipIndex();
 472
 473   Bool_t kPatRec = kFALSE;
 474
 475   Int_t CandidatePhotons = 0;
 476
 477   Float_t ShiftX = GetShiftX();
 478   Float_t ShiftY = GetShiftY();
 479
 480   Float_t* CandidatePhotonX = GetCandidatePhotonX();
 481   Float_t* CandidatePhotonY = GetCandidatePhotonY();
 482
 483   Int_t CandidatePhotonsNumber = GetCandidatePhotonsNumber();
 484
 485   if(fDebug) cout << " n " << CandidatePhotonsNumber << endl;
 486
 487   SetThetaCerenkov(999.);
 488   SetHoughPhotons(0);
 489   SetHoughPhotonsNorm(0);
 490   SetHoughRMS(999.);
 491
 492   for (Int_t j=0; j < CandidatePhotonsNumber; j++)
 493     {
 494
 495       SetPhotonIndex(j);
 496
 497       SetPhotonFlag(0);
 498       SetPhotonEta(-999.);
 499       SetPhotonWeight(0.);
 500
 501       if (j == MipIndex) continue;
 502
 503       //      h2_test1->Fill(CandidatePhotonX[j],CandidatePhotonY[j]);
 504
 505       if(CandidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */
 506
 507       Float_t Xtoentr = CandidatePhotonX[j] - ShiftX;
 508       Float_t Ytoentr = CandidatePhotonY[j] - ShiftY;
 509
 510       //      Float_t chargehit = fHits_charge[j];
 511       //      if(chargehit > 150) continue;
 512
 513       SetEntranceX(Xtoentr);
 514       SetEntranceY(Ytoentr);
 515
 516       FindPhiPoint();
 517
 518       Int_t PhotonStatus = PhotonInBand();
 519
 520       if(fDebug)
 521          {
 522             cout << " Photon n. " << j << " Status " << PhotonStatus << " accepted " << endl;
 523             cout << " CandidatePhotonX[j] " << CandidatePhotonX[j] << " CandidatePhotonY[j] " << CandidatePhotonY[j] << endl;
 524          }
 525
 526       if(PhotonStatus == 0) continue;
 527
 528       SetPhotonFlag(1);
 529
 530       FindThetaPhotonCerenkov();
 531
 532       Float_t ThetaPhotonCerenkov = GetThetaPhotonCerenkov();
 533
 534       if(fDebug) cout << " theta photon " << ThetaPhotonCerenkov << endl;
 535
 536       SetPhotonEta(ThetaPhotonCerenkov);
 537
 538       CandidatePhotons++;
 539
 540       // fill histograms
 541
 542       //      h2_test4->Fill(CandidatePhotonX[j],CandidatePhotonY[j]);
 543
 544       //      if(kDISPLAY) h1_photons->Fill(ThetaPhotonCerenkov);
 545
 546     }
 547
 548   if(CandidatePhotons >= 1) kPatRec = kTRUE;
 549
 550   if(!kPatRec) return;
 551     {
 552        SetThetaCerenkov(999.);
 553        SetHoughPhotons(0);
 554     }
 555   SetPhotonsNumber(CandidatePhotonsNumber);
 556
 557   HoughResponse();
 558
 559   NRings++;
 560
 561   FlagPhotons();
 562   Int_t NPhotonHough = GetHoughPhotons();
 563
 564   if(NPhotonHough < 1)
 565     {
 566       SetThetaCerenkov(999.);
 567       SetHoughPhotonsNorm(0.);
 568       return;
 569     }
 570
 571   if(kWEIGHT) FindWeightThetaCerenkov();
 572
 573   Float_t ThetaCerenkov = GetThetaCerenkov();
 574
 575   SetThetaOfRing(ThetaCerenkov);
 576   FindAreaAndPortionOfRing();
 577
 578   Float_t NPhotonHoughNorm = ((Float_t)NPhotonHough)/GetPortionOfRing();
 579   SetHoughPhotonsNorm(NPhotonHoughNorm);
 580
 581   // Calculate the area where the photon are accepted...
 582
 583   Float_t ThetaInternal = ThetaCerenkov - 0.5*fWindowWidth;
 584   SetThetaOfRing(ThetaInternal);
 585   FindAreaAndPortionOfRing();
 586   Float_t InternalArea = GetAreaOfRing();
 587
 588   Float_t ThetaExternal = ThetaCerenkov + 0.5*fWindowWidth;
 589   SetThetaOfRing(ThetaExternal);
 590   FindAreaAndPortionOfRing();
 591   Float_t ExternalArea = GetAreaOfRing();
 592
 593   Float_t HoughArea = ExternalArea - InternalArea;
 594
 595   SetHoughArea(HoughArea);
 596
 597   if(fDebug)
 598     {
 599       cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl;
 600       cout << " Rings found " << NRings << " with thetac " << ThetaCerenkov << endl;
 601
 602       h1_hough->Fill(ThetaCerenkov,1.);
 603
 604       cout << " Nphotons " << GetPhotonsNumber()
 605            << " Hough    " << NPhotonHough
 606            << " norm     " << NPhotonHoughNorm << endl;
 607
 608       cout << " In PatRec:p " << pmod << " theta " << TrackTheta << " phi " << TrackPhi << endl;
 609       cout << " ------------------------------------- " << endl;
 610     }
 611
 612   Int_t NPhotons = GetPhotonsNumber();
 613
 614   Float_t xmean = 0.;
 615   Float_t x2mean = 0.;
 616   Int_t nev = 0;
 617
 618   for (Int_t j=0; j < NPhotons;j++)
 619     {
 620       SetPhotonIndex(j);
 621
 622       Float_t eta = GetPhotonEta();
 623
 624       if(eta != -999.)
 625         {
 626           if(GetPhotonFlag() == 2)
 627             {
 628
 629               if(pmod>2.5&&ThetaCerenkov>0.65) photris->Fill(eta);
 630
 631               xmean += eta;
 632               x2mean += eta*eta;
 633               nev++;
 634             }
 635         }
 636     }
 637
 638   if(nev > 0)
 639     {
 640       xmean /=(Float_t)nev;
 641       x2mean /=(Float_t)nev;
 642     } else {
 643       xmean = 0.;
 644       x2mean = 0.;
 645     }
 646
 647   Float_t RMS = sqrt(x2mean - xmean*xmean);
 648
 649   SetHoughRMS(RMS);
 650
 651   if(fDebug) cout << " RMS " << RMS << endl;
 652
 653 }
 654
 655 void AliRICHRecon::FindEmissionPoint()
 656 {
 657
 658 // Find emission point
 659
 660   Float_t absorbtionLenght=7.83*RadiatorWidth; //absorption length in the freon (cm)
 661   // 7.83 = -1/ln(T0) where
 662   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
 663   Float_t photonLenght, photonLenghtMin, photonLenghtMax;
 664
 665   photonLenght=exp(-RadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
 666   photonLenghtMin=RadiatorWidth*photonLenght/(1.-photonLenght);
 667   photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
 668   Float_t EmissionPoint = RadiatorWidth + photonLenghtMin - photonLenghtMax;
 669
 670   SetEmissionPoint(EmissionPoint);
 671 }
 672
 673
 674 Int_t AliRICHRecon::PhotonInBand()
 675 {
 676
 677   //  Float_t MassOfParticle;
 678   Float_t beta;
 679   Float_t nfreon;
 680
 681   Float_t thetacer;
 682
 683   Float_t Xtoentr = GetEntranceX();
 684   Float_t Ytoentr = GetEntranceY();
 685
 686   Float_t InnerRadius;
 687   Float_t OuterRadius;
 688
 689   Float_t phpad = GetPhiPoint();
 690
 691   //  Float_t pmod = GetTrackMomentum();
 692   //  Float_t TrackTheta = GetTrackTheta();
 693   //  Float_t TrackPhi = GetTrackPhi();
 694
 695   // inner radius //
 696   SetPhotonEnergy(5.6);
 697   SetEmissionPoint(RadiatorWidth -0.0001);
 698   SetMassHypotesis(0.93828);
 699
 700   SetBetaOfParticle();
 701   SetFreonRefractiveIndex();
 702
 703   beta   = GetBetaOfParticle();
 704   nfreon = GetFreonRefractiveIndex();
 705
 706   thetacer = Cerenkovangle(nfreon,beta);
 707
 708   thetacer = 0.;
 709
 710   if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl;
 711
 712   FindThetaAtQuartz(thetacer);
 713
 714   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 715     {
 716       InnerRadius = -999.;
 717       SetXInnerRing(-999.);
 718       SetYInnerRing(-999.);
 719       SetRadiusInnerRing(-999.);
 720     }
 721   else
 722     {
 723       SetThetaPhotonInDRS(GetThetaAtQuartz());
 724       SetPhiPhotonInDRS(phpad);
 725
 726       InnerRadius = FromEmissionToCathode();
 727        if(InnerRadius == 999.) InnerRadius = -999.;
 728
 729       SetXInnerRing(GetXPointOnCathode());
 730       SetYInnerRing(GetYPointOnCathode());
 731       SetRadiusInnerRing(InnerRadius);
 732     }
 733
 734   // outer radius //
 735   SetPhotonEnergy(7.7);
 736   SetEmissionPoint(0.);
 737 //  SetMassHypotesis(0.139567);
 738   SetMassHypotesis(0.);
 739
 740   SetBetaOfParticle();
 741   SetFreonRefractiveIndex();
 742
 743   beta   = GetBetaOfParticle();
 744   nfreon = GetFreonRefractiveIndex();
 745
 746   thetacer = Cerenkovangle(nfreon,beta);
 747
 748   //  thetacer = 0.75;
 749
 750   if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl;
 751
 752   FindThetaAtQuartz(thetacer);
 753
 754   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 755     {
 756       OuterRadius = 999.;
 757       SetXOuterRing(999.);
 758       SetYOuterRing(999.);
 759       SetRadiusOuterRing(999.);
 760     }
 761   else
 762     {
 763       SetThetaPhotonInDRS(GetThetaAtQuartz());
 764       SetPhiPhotonInDRS(phpad);
 765
 766       OuterRadius = FromEmissionToCathode();
 767 //      cout << " OuterRadius " << OuterRadius << endl;
 768       SetXOuterRing(GetXPointOnCathode());
 769       SetYOuterRing(GetYPointOnCathode());
 770       SetRadiusOuterRing(OuterRadius);
 771     }
 772
 773   Float_t padradius = sqrt(TMath::Power(Xtoentr,2)+TMath::Power(Ytoentr,2));
 774
 775   if(fDebug) printf(" rmin %f r %f rmax %f \n",InnerRadius,padradius,OuterRadius);
 776
 777   if(padradius>=InnerRadius && padradius<=OuterRadius) return 1;
 778   return 0;
 779 }
 780
 781 void AliRICHRecon::FindThetaAtQuartz(Float_t ThetaCerenkov)
 782 {
 783
 784   if(ThetaCerenkov == 999.)
 785     {
 786       SetThetaAtQuartz(999.);
 787       return;
 788     }
 789
 790   Float_t ThetaAtQuartz = 999.;
 791
 792   Float_t TrackTheta = GetTrackTheta();
 793
 794   if(TrackTheta == 0) {
 795
 796     if(fDebug) cout << " Theta sol unique " << ThetaCerenkov << endl;
 797
 798     ThetaAtQuartz = ThetaCerenkov;
 799     SetThetaAtQuartz(ThetaAtQuartz);
 800     return;
 801   }
 802
 803   Float_t TrackPhi   = GetTrackPhi();
 804   Float_t PhiPoint = GetPhiPoint();
 805
 806   Double_t den = TMath::Sin((Double_t)TrackTheta)
 807     *TMath::Cos((Double_t)TrackPhi)
 808     *TMath::Cos((Double_t)PhiPoint) +
 809     TMath::Sin((Double_t)TrackTheta)
 810     *TMath::Sin((Double_t)TrackPhi)
 811     *TMath::Sin((Double_t)PhiPoint);
 812   Double_t b = TMath::Cos((Double_t)TrackTheta)/den;
 813   Double_t c = -TMath::Cos((Double_t)ThetaCerenkov)/den;
 814
 815   Double_t UnderSqrt = 1 + b*b - c*c;
 816
 817   if(fDebug)
 818     {
 819       cout << " TrackTheta    " << TrackTheta    << endl;
 820       cout << " TrackPhi      " << TrackPhi      << endl;
 821       cout << " PhiPoint      " << PhiPoint      << endl;
 822       cout << " ThetaCerenkov " << ThetaCerenkov << endl;
 823       cout << " den b c " << den << " b " << b << " c " << c << endl;
 824     }
 825
 826   if(UnderSqrt < 0) {
 827     if(fDebug) cout << " sqrt negative !!!!" << UnderSqrt << endl;
 828     SetThetaAtQuartz(999.);
 829     return;
 830   }
 831
 832   Double_t sol1 = (1+TMath::Sqrt(UnderSqrt))/(b-c);
 833   Double_t sol2 = (1-TMath::Sqrt(UnderSqrt))/(b-c);
 834
 835   Double_t ThetaSol1 = 2*TMath::ATan(sol1);
 836   Double_t ThetaSol2 = 2*TMath::ATan(sol2);
 837
 838   if(fDebug) cout << " Theta sol 1 " << ThetaSol1
 839                   << " Theta sol 2 " << ThetaSol2 << endl;
 840
 841   if(ThetaSol1>0 && ThetaSol1 < pi) ThetaAtQuartz = (Float_t)ThetaSol1;
 842   if(ThetaSol2>0 && ThetaSol2 < pi) ThetaAtQuartz = (Float_t)ThetaSol2;
 843
 844   SetThetaAtQuartz(ThetaAtQuartz);
 845 }
 846
 847 void AliRICHRecon::FindThetaPhotonCerenkov()
 848 {
 849
 850   Float_t ThetaCerMin = 0.;
 851   Float_t ThetaCerMax = 0.75;
 852   Float_t ThetaCerMean;
 853
 854   Float_t RadiusMin, RadiusMax, RadiusMean;
 855   Int_t nIteration = 0;
 856
 857   const Float_t Tollerance = 0.05;
 858
 859   //  Float_t pmod = GetTrackMomentum();
 860   //  Float_t TrackTheta = GetTrackTheta();
 861   //  Float_t TrackPhi = GetTrackPhi();
 862
 863   Float_t PhiPoint = GetPhiPoint();
 864
 865   SetPhotonEnergy(6.85);
 866   SetEmissionPoint(RadiatorWidth/2);
 867
 868   Float_t XPoint = GetEntranceX();
 869   Float_t YPoint = GetEntranceY();
 870   Float_t DistPoint = sqrt(XPoint*XPoint + YPoint*YPoint);
 871
 872   if(fDebug) cout << " DistPoint " << DistPoint << endl;
 873
 874   // Star minimization...
 875
 876   // First value...
 877
 878   FindThetaAtQuartz(ThetaCerMin);
 879
 880   if(GetThetaAtQuartz() == 999.)
 881     {
 882       RadiusMin = -999.;
 883     }
 884   else
 885     {
 886       SetThetaPhotonInDRS(GetThetaAtQuartz());
 887       SetPhiPhotonInDRS(PhiPoint);
 888
 889       RadiusMin = FromEmissionToCathode();
 890     }
 891
 892   // Second value...
 893
 894   FindThetaAtQuartz(ThetaCerMax);
 895   if(GetThetaAtQuartz() == 999.)
 896     {
 897       RadiusMax = 999.;
 898     }
 899   else
 900     {
 901       SetThetaPhotonInDRS(GetThetaAtQuartz());
 902       SetPhiPhotonInDRS(PhiPoint);
 903
 904       RadiusMax = FromEmissionToCathode();
 905     }
 906   // Mean value...
 907
 908   ThetaCerMean = (ThetaCerMax + ThetaCerMin)/2;
 909
 910   FindThetaAtQuartz(ThetaCerMean);
 911   if(GetThetaAtQuartz() == 999.)
 912     {
 913       RadiusMean = 999.;
 914     }
 915   else
 916     {
 917       SetThetaPhotonInDRS(GetThetaAtQuartz());
 918       SetPhiPhotonInDRS(PhiPoint);
 919
 920       RadiusMean = FromEmissionToCathode();
 921     }
 922
 923   if(fDebug) cout << " r1 " << RadiusMin << " rmean "
 924                   << RadiusMean << " r2 " << RadiusMax << endl;
 925
 926   while (TMath::Abs(RadiusMean-DistPoint) > Tollerance)
 927     {
 928
 929       if((RadiusMin-DistPoint)*(RadiusMean-DistPoint) < 0) ThetaCerMax = ThetaCerMean;
 930       if((RadiusMin-DistPoint)*(RadiusMean-DistPoint) > 0) {
 931
 932         ThetaCerMin = ThetaCerMean;
 933
 934         FindThetaAtQuartz(ThetaCerMin);
 935         SetThetaPhotonInDRS(GetThetaAtQuartz());
 936         SetPhiPhotonInDRS(PhiPoint);
 937
 938         RadiusMin =FromEmissionToCathode();
 939       }
 940
 941       ThetaCerMean = (ThetaCerMax + ThetaCerMin)/2;
 942
 943       FindThetaAtQuartz(ThetaCerMean);
 944       SetThetaPhotonInDRS(GetThetaAtQuartz());
 945       SetPhiPhotonInDRS(PhiPoint);
 946
 947       RadiusMean = FromEmissionToCathode();
 948
 949       nIteration++;
 950       if(nIteration>=50) {
 951         if(fDebug) printf(" max iterations in FindPhotonCerenkov\n");
 952         SetThetaPhotonCerenkov(999.);
 953         return;
 954       }
 955     }
 956
 957   SetThetaPhotonCerenkov(ThetaCerMean);
 958
 959 }
 960
 961 void AliRICHRecon::FindAreaAndPortionOfRing()
 962 {
 963
 964   Float_t XPoint[NPointsOfRing], YPoint[NPointsOfRing];
 965
 966   //  Float_t Xtoentr = GetEntranceX();
 967   //  Float_t Ytoentr = GetEntranceY();
 968   Float_t ShiftX = GetShiftX();
 969   Float_t ShiftY = GetShiftY();
 970
 971   Float_t XEmiss = GetXCoordOfEmission();
 972   Float_t YEmiss = GetYCoordOfEmission();
 973
 974   Float_t x0 = XEmiss + ShiftX;
 975   Float_t y0 = YEmiss + ShiftY;
 976
 977   //  Float_t pmod = GetTrackMomentum();
 978   //  Float_t TrackTheta = GetTrackTheta();
 979   //  Float_t TrackPhi = GetTrackPhi();
 980
 981   SetPhotonEnergy(6.85);
 982   SetFreonRefractiveIndex();
 983
 984   SetEmissionPoint(RadiatorWidth/2.);
 985
 986   Float_t Theta = GetThetaOfRing();
 987
 988   Int_t nPoints = 0;
 989   Int_t NPsiAccepted = 0;
 990   Int_t NPsiTotal = 0;
 991
 992   for(Int_t i=0;i<NPointsOfRing-1;i++)
 993     {
 994
 995       Float_t Psi = 2*TMath::Pi()*i/NPointsOfRing;
 996
 997       SetThetaPhotonInTRS(Theta);
 998       SetPhiPhotonInTRS(Psi);
 999       FindPhotonAnglesInDRS();
1000
1001       Float_t Radius = FromEmissionToCathode();
1002       if (Radius == 999.) continue;
1003
1004       NPsiTotal++;
1005
1006       Float_t XPointRing = GetXPointOnCathode() + ShiftX;
1007       Float_t YPointRing = GetYPointOnCathode() + ShiftY;
1008
1009       SetDetectorWhereX(XPointRing);
1010       SetDetectorWhereY(YPointRing);
1011
1012       Int_t Zone = CheckDetectorAcceptance();
1013
1014 //      cout << " XPointing " << XPointRing << " YPointing " << YPointRing << " Zone " << Zone << endl;
1015 //      cout << " ShiftX " << ShiftX << " ShiftY " << ShiftY << endl;
1016 //      cout << " GetXPointOnCathode() " << GetXPointOnCathode() << endl;
1017 //      cout << " GetYPointOnCathode() " << GetYPointOnCathode() << endl;
1018
1019       if (Zone != 0)
1020         {
1021           FindIntersectionWithDetector();
1022           XPoint[nPoints] = GetIntersectionX();
1023           YPoint[nPoints] = GetIntersectionY();
1024         }
1025       else
1026         {
1027           XPoint[nPoints] = XPointRing;
1028           YPoint[nPoints] = YPointRing;
1029           NPsiAccepted++;
1030         }
1031
1032       nPoints++;
1033
1034     }
1035
1036   XPoint[nPoints] = XPoint[0];
1037   YPoint[nPoints] = YPoint[0];
1038
1039   // find area...
1040
1041   Float_t Area = 0;
1042
1043   for (Int_t i = 0; i < nPoints; i++)
1044     {
1045       Area += TMath::Abs((XPoint[i]-x0)*(YPoint[i+1]-y0) - (XPoint[i+1]-x0)*(YPoint[i]-y0));
1046     }
1047
1048   Area *= 0.5;
1049
1050   Float_t PortionOfRing = ((Float_t)NPsiAccepted)/((Float_t)(NPsiTotal));
1051
1052   //  cout << " Area " << Area << " Portion of ring " << PortionOfRing << endl;
1053
1054   SetAreaOfRing(Area);
1055   SetPortionOfRing(PortionOfRing);
1056 }
1057
1058 void AliRICHRecon::FindIntersectionWithDetector()
1059 {
1060
1061   Float_t XIntersect, YIntersect;
1062   Float_t x1, x2, y1, y2;
1063
1064   Float_t ShiftX = GetShiftX();
1065   Float_t ShiftY = GetShiftY();
1066
1067   Float_t XPoint = GetXPointOnCathode() + ShiftX;
1068   Float_t YPoint = GetYPointOnCathode() + ShiftY;
1069
1070   Float_t XEmiss = GetXCoordOfEmission();
1071   Float_t YEmiss = GetYCoordOfEmission();
1072
1073   Float_t Phi = GetPhiPhotonInDRS();
1074   Float_t m = tan(Phi);
1075
1076   Float_t x0 = XEmiss + ShiftX;
1077   Float_t y0 = YEmiss + ShiftY;
1078
1079   if(XPoint > x0)
1080     {
1081       x1 = x0;
1082       x2 = XPoint;
1083     }
1084   else
1085     {
1086       x2 = x0;
1087       x1 = XPoint;
1088     }
1089   if(YPoint > y0)
1090     {
1091       y1 = y0;
1092       y2 = YPoint;
1093     }
1094   else
1095     {
1096       y2 = y0;
1097       y1 = YPoint;
1098     }
1099   //
1100   XIntersect = Xmax;
1101   YIntersect = m*(XIntersect - x0) + y0;
1102   if (YIntersect >= Ymin && YIntersect <= Ymax && XIntersect >= x1 && XIntersect <= x2)
1103     {
1104       SetIntersectionX(XIntersect);
1105       SetIntersectionY(YIntersect);
1106       return;
1107     }
1108   //
1109   XIntersect = Xmin;
1110   YIntersect = m*(XIntersect - x0) + y0;
1111   if (YIntersect >= Ymin && YIntersect <= Ymax && XIntersect >= x1 && XIntersect <= x2)
1112     {
1113       SetIntersectionX(XIntersect);
1114       SetIntersectionY(YIntersect);
1115       return;
1116     }
1117   //
1118   YIntersect = Ymax;
1119   XIntersect = (YIntersect - y0)/m + x0;
1120   if (XIntersect >= Xmin && XIntersect <= Xmax && YIntersect >= y1 && YIntersect <= y2)
1121     {
1122       SetIntersectionX(XIntersect);
1123       SetIntersectionY(YIntersect);
1124       return;
1125     }
1126   //
1127   YIntersect = Ymin;
1128   XIntersect = (YIntersect - y0)/m + x0;
1129   if (XIntersect >= Xmin && XIntersect <= Xmax && YIntersect >= y1 && YIntersect <= y2)
1130     {
1131       SetIntersectionX(XIntersect);
1132       SetIntersectionY(YIntersect);
1133       return;
1134     }
1135
1136   cout << " sono fuori!!!!!!" << endl;
1137 //  cout << " x1 " << x1 << " x2 " << x2 << endl;
1138 //  cout << " y1 " << y1 << " y2 " << y2 << endl;
1139 //  cout << " Xmin " << Xmin << " Xmax " << Xmax << endl;
1140 //  cout << " Ymin " << Ymin << " Ymax " << Ymax << endl;
1141
1142 }
1143
1144 Int_t AliRICHRecon::CheckDetectorAcceptance()
1145 {
1146
1147   // crosses X -2.6 2.6 cm
1148   // crosses Y -1 1 cm
1149
1150   Float_t Xcoord = GetDetectorWhereX();
1151   Float_t Ycoord = GetDetectorWhereY();
1152
1153 //  cout << " Xcoord " << Xcoord << " Ycoord " << Ycoord << endl;
1154   if(Xcoord > Xmax)
1155     {
1156       if(Ycoord > Ymax) return 2;
1157       if(Ycoord > Ymin && Ycoord < Ymax) return 3;
1158       if(Ycoord < Ymin) return 4;
1159     }
1160   if(Xcoord < Xmin)
1161     {
1162       if(Ycoord > Ymax) return 8;
1163       if(Ycoord > Ymin && Ycoord < Ymax) return 7;
1164       if(Ycoord < Ymin) return 6;
1165     }
1166   if(Xcoord > Xmin && Xcoord < Xmax)
1167     {
1168       if(Ycoord > Ymax) return 1;
1169       if(Ycoord > Ymin && Ycoord < Ymax) return 0;
1170       if(Ycoord < Ymin) return 5;
1171     }
1172   return 999;
1173 }
1174
1175 void AliRICHRecon::DrawRing()
1176 {
1177
1178   //  Float_t xGraph[1000],yGraph[1000];
1179
1180   Float_t type;
1181   //  Float_t MassOfParticle;
1182   Float_t beta;
1183   Float_t nfreon;
1184
1185   Float_t ThetaCerenkov;
1186
1187   //  Float_t Xtoentr = GetEntranceX();
1188   //  Float_t Ytoentr = GetEntranceY();
1189
1190   //  Float_t pmod = GetTrackMomentum();
1191   //  Float_t TrackTheta = GetTrackTheta();
1192   //  Float_t TrackPhi = GetTrackPhi();
1193
1194   SetPhotonEnergy(6.85);
1195   SetFreonRefractiveIndex();
1196
1197   SetEmissionPoint(RadiatorWidth/2.);
1198
1199   type = 1;
1200
1201   if(type == 1)
1202     {
1203       SetMassHypotesis(0.139567);
1204       SetBetaOfParticle();
1205
1206       beta   = GetBetaOfParticle();
1207
1208     }
1209   else if(type == 2)
1210     {
1211       ThetaCerenkov = GetThetaCerenkov();
1212       FindBetaFromTheta(ThetaCerenkov);
1213     }
1214
1215   nfreon = GetFreonRefractiveIndex();
1216
1217   Float_t thetacer = Cerenkovangle(nfreon,beta);
1218
1219   if(fDebug) cout << " TetaCer in DrawRing " << thetacer << endl;
1220
1221   Int_t nPoints = 100;
1222
1223   Int_t nPointsToDraw = 0;
1224   for(Int_t i=0;i<nPoints;i++)
1225     {
1226       Float_t phpad = 2*TMath::Pi()*i/nPoints;
1227       SetThetaPhotonInTRS(thetacer);
1228       SetPhiPhotonInTRS(phpad);
1229       FindPhotonAnglesInDRS();
1230       Float_t Radius = FromEmissionToCathode();
1231       if (Radius == 999.) continue;
1232       xGraph[nPointsToDraw] = GetXPointOnCathode() + GetShiftX();
1233       yGraph[nPointsToDraw] = GetYPointOnCathode() + GetShiftY();
1234       //      cout << " get shift X " << GetShiftX() << endl;
1235       //      cout << " get shift Y " << GetShiftY() << endl;
1236       nPointsToDraw++;
1237     }
1238
1239
1240   if(fDebug) cout << " Drawing the Ring... with " << nPointsToDraw << " points " << endl;
1241
1242   //  pol = new TPolyLine(nPointsToDraw,xGraph,yGraph);
1243   //  pol->Draw("same");
1244   gra = new TGraph(nPointsToDraw,xGraph,yGraph);
1245   gra->Draw("AC");
1246   StarCanvas->Update();
1247
1248 }
1249
1250 Float_t AliRICHRecon::PhotonPositionOnCathode()
1251 {
1252   //  Float_t MassOfParticle;
1253   Float_t beta;
1254   Float_t nfreon;
1255
1256   //  Float_t pmod = GetTrackMomentum();
1257   //  Float_t TrackTheta = GetTrackTheta();
1258   //  Float_t TrackPhi = GetTrackPhi();
1259
1260   //  Float_t phpad = GetPhiPoint();
1261
1262   SetPhotonEnergy(6.85);
1263   SetEmissionPoint(RadiatorWidth/2.);
1264   SetMassHypotesis(0.139567);
1265
1266   SetBetaOfParticle();
1267   SetFreonRefractiveIndex();
1268
1269   beta   = GetBetaOfParticle();
1270   nfreon = GetFreonRefractiveIndex();
1271
1272   //  Float_t thetacer = Cerenkovangle(nfreon,beta);
1273
1274   //  cout << " FromEmissionToCathode: thetacer " << thetacer << " phpad " << phpad << endl;
1275
1276   Float_t Radius = FromEmissionToCathode();
1277   if (Radius == 999.) return 999.;
1278
1279   //  Float_t Xphoton = GetXPointOnCathode();
1280   //  Float_t Yphoton = GetYPointOnCathode();
1281   //  cout << " PhotonPositionOnCathode: Xphoton " << Xphoton << " Yphoton " << Yphoton <<
1282   //  " Radius for photon " << Radius << endl;
1283   return 0;
1284 }
1285
1286 void AliRICHRecon::FindPhotonAnglesInDRS()
1287 {
1288   // Setup the rotation matrix of the track...
1289
1290   TRotation Mtheta;
1291   TRotation Mphi;
1292   TRotation Minv;
1293   TRotation Mrot;
1294
1295   Float_t TrackTheta = GetTrackTheta();
1296   Float_t TrackPhi = GetTrackPhi();
1297
1298   Mtheta.RotateY(TrackTheta);
1299   Mphi.RotateZ(TrackPhi);
1300
1301   Mrot = Mphi * Mtheta;
1302   //  Minv = Mrot.Inverse();
1303
1304   TVector3 PhotonInRadiator(1,1,1);
1305
1306   Float_t ThetaCerenkov = GetThetaPhotonInTRS();
1307   Float_t PhiCerenkov   = GetPhiPhotonInTRS();
1308
1309   PhotonInRadiator.SetTheta(ThetaCerenkov);
1310   PhotonInRadiator.SetPhi(PhiCerenkov);
1311   PhotonInRadiator = Mrot * PhotonInRadiator;
1312   Float_t Theta = PhotonInRadiator.Theta();
1313   Float_t Phi = PhotonInRadiator.Phi();
1314   SetThetaPhotonInDRS(Theta);
1315   SetPhiPhotonInDRS(Phi);
1316
1317 }
1318
1319 Float_t AliRICHRecon::FromEmissionToCathode()
1320 {
1321
1322   Float_t nfreon, nquartz, ngas;
1323
1324   SetFreonRefractiveIndex();
1325   SetQuartzRefractiveIndex();
1326   SetGasRefractiveIndex();
1327
1328   nfreon  = GetFreonRefractiveIndex();
1329   nquartz = GetQuartzRefractiveIndex();
1330   ngas    = GetGasRefractiveIndex();
1331
1332   Float_t TrackTheta = GetTrackTheta();
1333   Float_t TrackPhi = GetTrackPhi();
1334   Float_t LengthOfEmissionPoint = GetEmissionPoint();
1335
1336   Float_t Theta = GetThetaPhotonInDRS();
1337   Float_t Phi   = GetPhiPhotonInDRS();
1338
1339 //   cout << " Theta " << Theta << " Phi " << Phi << endl;
1340
1341   Float_t xEmiss = LengthOfEmissionPoint*tan(TrackTheta)*cos(TrackPhi);
1342   Float_t yEmiss = LengthOfEmissionPoint*tan(TrackTheta)*sin(TrackPhi);
1343
1344   SetXCoordOfEmission(xEmiss);
1345   SetYCoordOfEmission(yEmiss);
1346
1347   Float_t thetaquar = SnellAngle(nfreon, nquartz, Theta);
1348
1349   if(thetaquar == 999.)
1350     {
1351       SetXPointOnCathode(999.);
1352       SetYPointOnCathode(999.);
1353       return thetaquar;
1354     }
1355
1356   Float_t thetagap  = SnellAngle( nquartz, ngas, thetaquar);
1357
1358   if(thetagap == 999.)
1359     {
1360       SetXPointOnCathode(999.);
1361       SetYPointOnCathode(999.);
1362       return thetagap;
1363     }
1364
1365   Float_t xw = (RadiatorWidth - LengthOfEmissionPoint)*cos(Phi)*tan(Theta);
1366   Float_t xq = QuartzWidth*cos(Phi)*tan(thetaquar);
1367   Float_t xg = GapWidth*cos(Phi)*tan(thetagap);
1368   Float_t yw = (RadiatorWidth - LengthOfEmissionPoint)*sin(Phi)*tan(Theta);
1369   Float_t yq = QuartzWidth*sin(Phi)*tan(thetaquar);
1370   Float_t yg = GapWidth*sin(Phi)*tan(thetagap);
1371
1372 //  Float_t xtot = x1 + xw + xq + xg;
1373 //  Float_t ytot = y1 + yw + yq + yg;
1374
1375   Float_t xtot = xEmiss + xw + xq + xg;
1376   Float_t ytot = yEmiss + yw + yq + yg;
1377
1378   SetXPointOnCathode(xtot);
1379   SetYPointOnCathode(ytot);
1380
1381 //  cout << " xtot " << xtot << " ytot " << ytot << endl;
1382
1383   Float_t DistanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
1384                                     +TMath::Power(fPhotonLimitY,2));
1385
1386   return DistanceFromEntrance;
1387
1388 }
1389
1390
1391 void AliRICHRecon::FindPhiPoint()
1392 {
1393
1394   Float_t Xtoentr = GetEntranceX();
1395   Float_t Ytoentr = GetEntranceY();
1396
1397   Float_t TrackTheta = GetTrackTheta();
1398   Float_t TrackPhi = GetTrackPhi();
1399
1400   Float_t EmissionPoint = GetEmissionPoint();
1401
1402   Float_t argY = Ytoentr - EmissionPoint*tan(TrackTheta)*sin(TrackPhi);
1403   Float_t argX = Xtoentr - EmissionPoint*tan(TrackTheta)*cos(TrackPhi);
1404   Float_t phipad = atan2(argY,argX);
1405
1406   SetPhiPoint(phipad);
1407
1408 }
1409
1410 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
1411 {
1412
1413 // Compute the cerenkov angle
1414
1415   Float_t thetacer;
1416
1417   if((n*beta)<1.) {
1418     thetacer = 999.;
1419     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
1420     return thetacer;
1421   }
1422
1423   thetacer = acos (1./(n*beta));
1424   return thetacer;
1425 }
1426
1427 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
1428 {
1429
1430 // Compute the Snell angle
1431
1432   Float_t sinrefractangle;
1433   Float_t refractangle;
1434
1435   sinrefractangle = (n1/n2)*sin(theta1);
1436
1437   if(sinrefractangle>1.) {
1438     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
1439     refractangle = 999.;
1440     return refractangle;
1441   }
1442
1443   refractangle = asin(sinrefractangle);
1444   return refractangle;
1445 }
1446
1447
1448 void AliRICHRecon::HoughResponse()
1449
1450 {
1451
1452 // Implement Hough response pat. rec. method
1453
1454   Float_t *HCSspace;
1455
1456   int           bin=0;
1457   int           bin1=0;
1458   int           bin2=0;
1459   int           i, j, k, nCorrBand;
1460   float         hcs[750],hcsw[750];
1461   float         angle, weight;
1462   float         lowerlimit,upperlimit;
1463
1464   float         etaPeak[100];
1465
1466   int           nBin;
1467
1468   float etaPeakPos  = -1;
1469
1470   Int_t   etaPeakCount = -1;
1471
1472   Float_t ThetaCerenkov = 0.;
1473
1474   nBin = (int)(0.5+fThetaMax/(fDTheta));
1475   nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta));
1476
1477   memset ((void *)hcs, 0, fThetaBin*sizeof(float));
1478   memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
1479
1480   Int_t NPhotons = GetPhotonsNumber();
1481
1482   Int_t WeightFlag = 0;
1483
1484   for (k=0; k< NPhotons; k++) {
1485
1486     SetPhotonIndex(k);
1487
1488     angle = GetPhotonEta();
1489
1490     if(angle == -999.) continue;
1491
1492     if (angle>=fThetaMin && angle<= fThetaMax)
1493
1494       {
1495
1496         bin = (int)(0.5+angle/(fDTheta));
1497
1498         bin1= bin-nCorrBand;
1499         bin2= bin+nCorrBand;
1500
1501         // calculate weights
1502
1503         if(kWEIGHT)
1504           {
1505             lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
1506             SetThetaOfRing(lowerlimit);
1507             FindAreaAndPortionOfRing();
1508             Float_t area1 = GetAreaOfRing();
1509
1510             upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
1511             SetThetaOfRing(upperlimit);
1512             FindAreaAndPortionOfRing();
1513             Float_t area2 = GetAreaOfRing();
1514
1515             //      cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
1516             Float_t diffarea = area2 - area1;
1517
1518             if(diffarea>0)
1519               {
1520                 weight = 1./(area2-area1);
1521               }
1522             else
1523               {
1524                 WeightFlag = 1;
1525                 weight = 1.;
1526               }
1527
1528             //      cout <<" low "<< lowerlimit << " up " << upperlimit <<
1529             //        " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
1530
1531           }
1532         else
1533           {
1534             weight = 1.;
1535           }
1536
1537         SetPhotonWeight(weight);
1538
1539         //      cout << "weight..." << weight << endl;
1540
1541         h1_photons1->Fill(angle);
1542         h1_photons2->Fill(angle,weight);
1543
1544         if (bin1<0)    bin1=0;
1545         if (bin2>nBin) bin2=nBin;
1546
1547         for (j=bin1; j<bin2; j++)
1548           {
1549             hcs[j] += 1;
1550             hcsw[j] += weight;
1551           }
1552       }
1553   }
1554
1555 //   if(kDISPLAY)
1556 //     {
1557 //       for(Int_t j=0;j<750;j++)
1558 //      {
1559 //        h1_hcs->Fill(((Float_t)j),hcs[j]);
1560 //        h1_hcsw->Fill(((Float_t)j),hcsw[j]);
1561 //      }
1562 //     }
1563
1564   if(WeightFlag == 0)
1565     {
1566       HCSspace = hcsw;
1567     }
1568   else
1569     {
1570       HCSspace = hcs;
1571       //      cout << " probems with weight...normal procedure adopted " << endl;
1572     }
1573
1574   HoughFiltering(HCSspace);
1575
1576   for (bin=0; bin <nBin; bin++) {
1577     angle = (bin+0.5) * (fDTheta);
1578     if (HCSspace[bin] && HCSspace[bin] > etaPeakPos) {
1579       etaPeakCount = 0;
1580       etaPeakPos = HCSspace[bin];
1581       etaPeak[0]=angle;
1582     }
1583     else {
1584       if (HCSspace[bin] == etaPeakPos) {
1585         etaPeak[++etaPeakCount] = angle;
1586       }
1587     }
1588   }
1589
1590   for (i=0; i<etaPeakCount+1; i++) {
1591     ThetaCerenkov += etaPeak[i];
1592   }
1593   if (etaPeakCount>=0) {
1594     ThetaCerenkov /= etaPeakCount+1;
1595     fThetaPeakPos = etaPeakPos;
1596   }
1597
1598   SetThetaCerenkov(ThetaCerenkov);
1599 }
1600
1601
1602 void AliRICHRecon::HoughFiltering(float hcs[])
1603 {
1604
1605 // hough filtering
1606
1607    float hcsFilt[750];
1608    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
1609    int nx, i, nxDx;
1610    int sizeHCS;
1611    int nBin;
1612
1613    nBin =  (int)(1+fThetaMax/fDTheta);
1614    sizeHCS = fThetaBin*sizeof(float);
1615
1616    memset ((void *)hcsFilt, 0, sizeHCS);
1617
1618    for (nx = 0; nx < nBin; nx++) {
1619       for (i = 0; i < 5; i++)   {
1620         nxDx = nx + (i-2);
1621         if (nxDx> -1 && nxDx<nBin)
1622              hcsFilt[nx] +=  hcs[nxDx] * k[i];
1623       }
1624    }
1625
1626    for (nx = 0; nx < nBin; nx++) {
1627      hcs[nx] = hcsFilt[nx];
1628    }
1629 }
1630
1631 void AliRICHRecon::FindWeightThetaCerenkov()
1632 {
1633
1634   Float_t wei = 0.;
1635   Float_t WeightThetaCerenkov = 0.;
1636
1637   Int_t NPhotons = GetPhotonsNumber();
1638   for(Int_t i=0;i<NPhotons;i++)
1639     {
1640       SetPhotonIndex(i);
1641
1642       if(GetPhotonFlag() == 2)
1643         {
1644           Float_t PhotonEta = GetPhotonEta();
1645           Float_t PhotonWeight = GetPhotonWeight();
1646           WeightThetaCerenkov += PhotonEta*PhotonWeight;
1647           wei += PhotonWeight;
1648         }
1649     }
1650
1651   if(wei != 0.)
1652     {
1653       WeightThetaCerenkov /= wei;
1654     }
1655   else
1656     {
1657       WeightThetaCerenkov = 0.;
1658     }
1659
1660   SetThetaCerenkov(WeightThetaCerenkov);
1661
1662   cout << " thetac weighted -> " << WeightThetaCerenkov << endl;
1663 }
1664
1665
1666 void AliRICHRecon::FlagPhotons()
1667 {
1668
1669   Int_t NPhotonHough = 0;
1670
1671   Float_t ThetaCerenkov = GetThetaCerenkov();
1672   if(fDebug) cout << " fThetaCerenkov " << ThetaCerenkov << endl;
1673
1674   Float_t ThetaDist= ThetaCerenkov - fThetaMin;
1675   Int_t steps = (Int_t)(ThetaDist / fDTheta);
1676
1677   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
1678   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
1679   Float_t tavg = 0.5*(tmin+tmax);
1680
1681   tmin = tavg - 0.5*fWindowWidth;
1682   tmax = tavg + 0.5*fWindowWidth;
1683
1684   if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl;
1685   if(fDebug) cout << " thetac " << ThetaCerenkov << endl;
1686
1687   //  Int_t CandidatePhotonsNumber = GetCandidatePhotonsNumber();
1688
1689   Int_t NPhotons = GetPhotonsNumber();
1690
1691   //  for(Int_t i=0;i<CandidatePhotonsNumber;i++)
1692
1693   for(Int_t i=0;i<NPhotons;i++)
1694     {
1695       SetPhotonIndex(i);
1696
1697       Float_t PhotonEta = GetPhotonEta();
1698
1699       if(PhotonEta == -999.) continue;
1700
1701       if(PhotonEta >= tmin && PhotonEta <= tmax)
1702         {
1703           SetPhotonFlag(2);
1704           NPhotonHough++;
1705         }
1706     }
1707   SetHoughPhotons(NPhotonHough);
1708 }
1709
1710 void AliRICHRecon::DrawEvent(Int_t flag)
1711 {
1712
1713   flag=1; // dummy to be removed...
1714 /*
1715   Float_t xGraph[3000],yGraph[3000];
1716
1717   Float_t ThetaCerenkov;
1718
1719   // Display event...
1720
1721   gStyle->SetPalette(1,0);
1722
1723   if(flag == 0)
1724     {
1725
1726       //      Display = new TCanvas("Display","Star Display",0,0,1200,750);
1727
1728       Display->ToggleEventStatus();
1729       Display->Modified()
1730
1731       text = new TText(0,0,"");
1732       text->SetTextFont(61);
1733       text->SetTextSize(0.03);
1734       text->SetTextAlign(22);
1735
1736       Display->Resize();
1737
1738       h2_disp->Reset();
1739
1740       for(Int_t j=1;j<=nPixels;j++)
1741         {
1742           Float_t xpad = fPixels_localX[j-1];
1743           Float_t ypad = fPixels_localY[j-1];
1744           h2_disp->Fill(xpad,ypad,fPixels_charge[j-1]);
1745         }
1746
1747       h2_disp->SetMaximum(200);
1748       //      h2_disp->SetMaximum(1);
1749       h2_disp->SetStats(0);
1750       h2_disp->Draw("colz");
1751
1752       for(Int_t i=0; i<nRichPrimaries;i++)
1753
1754         {
1755
1756           TrackPoints = new TMarker(fRichPrimaries_localPadX[i],
1757                                     fRichPrimaries_localPadY[i],3);
1758
1759           TrackPoints->SetMarkerSize(1.5);
1760
1761           Float_t pmod = sqrt(fRichPrimaries_localPadPx[i] * fRichPrimaries_localPadPx[i] +
1762                               fRichPrimaries_localPadPy[i] * fRichPrimaries_localPadPy[i] +
1763                               fRichPrimaries_localPadPz[i] * fRichPrimaries_localPadPz[i]);
1764
1765           if(pmod < 1) TrackPoints->SetMarkerColor(kBlue);
1766           if(pmod > 1 && pmod < 2) TrackPoints->SetMarkerColor(kGreen);
1767           if(pmod > 2) TrackPoints->SetMarkerColor(kRed);
1768
1769           TrackPoints->Draw();
1770
1771           line = new TLine(-0.13,-42.,-0.13,42.);
1772           line->Draw();
1773           line = new TLine(0.13,-42.,0.13,42.);
1774           line->Draw();
1775           line = new TLine(-64.,-0.13,64.,-0.13);
1776           line->Draw();
1777           line = new TLine(-64.,0.13,64.,0.13);
1778           line->Draw();
1779
1780         }
1781
1782       return;
1783
1784     }
1785
1786   //
1787
1788   // Draw rings...
1789
1790   //
1791
1792   //  Float_t Xtoentr = GetEntranceX();
1793   //  Float_t Ytoentr = GetEntranceY();
1794
1795   //  Float_t pmod = GetTrackMomentum();
1796   //  Float_t TrackTheta = GetTrackTheta();
1797   //  Float_t TrackPhi = GetTrackPhi();
1798
1799   SetPhotonEnergy(6.85);
1800   SetFreonRefractiveIndex();
1801
1802   SetEmissionPoint(RadiatorWidth/2.);
1803
1804   ThetaCerenkov = GetThetaCerenkov();
1805
1806   if (ThetaCerenkov == 999.) return;
1807
1808   Int_t nPointsToDraw = 0;
1809
1810   for(Int_t i=0;i<99;i++)
1811     {
1812       Float_t phpad = 2*TMath::Pi()*i/99;
1813       SetThetaPhotonInTRS(ThetaCerenkov);
1814       SetPhiPhotonInTRS(phpad);
1815       FindPhotonAnglesInDRS();
1816       Float_t Radius = FromEmissionToCathode();
1817
1818       if (Radius == 999.) continue;
1819
1820       Float_t ShiftX = GetShiftX();
1821       Float_t ShiftY = GetShiftY();
1822
1823       Float_t XPointRing = GetXPointOnCathode() + ShiftX;
1824       Float_t YPointRing = GetYPointOnCathode() + ShiftY;
1825
1826       SetDetectorWhereX(XPointRing);
1827       SetDetectorWhereY(YPointRing);
1828
1829       Int_t Zone = CheckDetectorAcceptance();
1830
1831       if (Zone != 0)
1832         {
1833           FindIntersectionWithDetector();
1834           xGraph[nPointsToDraw] = GetIntersectionX();
1835           yGraph[nPointsToDraw] = GetIntersectionY();
1836           nPointsToDraw++;
1837         }
1838       else
1839         {
1840           xGraph[nPointsToDraw] = GetXPointOnCathode() + GetShiftX();
1841           yGraph[nPointsToDraw] = GetYPointOnCathode() + GetShiftY();
1842           nPointsToDraw++;
1843         }
1844     }
1845
1846   xGraph[nPointsToDraw] = xGraph[0];
1847   yGraph[nPointsToDraw] = yGraph[0];
1848
1849   poll = new TPolyLine(nPointsToDraw+1,xGraph,yGraph);
1850   poll->SetLineColor(2);
1851   poll->SetLineWidth(3);
1852
1853   Display->Update();
1854
1855   //  waiting();
1856   poll->Draw();
1857
1858   for(Int_t j=0;j<nHits;j++)
1859     {
1860
1861       Float_t xhit = fHits_localX[j];
1862       Float_t yhit = fHits_localY[j];
1863
1864       SetPhotonIndex(j);
1865       Int_t FlagPhoton = GetPhotonFlag();
1866
1867 //       if(FlagPhoton >= 1)
1868 //      {
1869
1870 //        Photon = new TMarker(xhit,yhit,4);
1871 //        Photon->SetMarkerSize(1.5);
1872 //        Photon->Draw("same");
1873
1874 //      }
1875
1876
1877       if(FlagPhoton == 2)
1878         {
1879
1880           PhotonAcc = new TMarker(xhit,yhit,30);
1881           PhotonAcc->SetMarkerSize(1.5);
1882           PhotonAcc->SetMarkerColor(50);
1883           PhotonAcc->Draw("same");
1884
1885         }
1886     }
1887
1888   Display->Update();
1889
1890 //   waiting();
1891 //   h1_photons->Draw();
1892 //   Display->Update();
1893
1894 //   waiting();
1895 //   h1_photacc->Draw();
1896 //   Display->Update();
1897
1898 //   waiting();
1899
1900 //   Display->Update();
1901
1902 //   h1_photons->Reset();
1903 //   h1_photacc->Reset();
1904
1905 */
1906 }
1907
1908 Float_t  AliRICHRecon::FindMassOfParticle()
1909 {
1910
1911   Float_t pmod = GetTrackMomentum();
1912
1913   SetPhotonEnergy(6.85);
1914   SetFreonRefractiveIndex();
1915
1916   Float_t ThetaCerenkov = GetThetaCerenkov();
1917   FindBetaFromTheta(ThetaCerenkov);
1918
1919   Double_t beta = (Double_t)(GetBetaOfParticle());
1920   Double_t den = 1. - beta*beta;
1921   if(den<=0.) return 999.;
1922
1923   Double_t gamma = 1./TMath::Sqrt(den);
1924
1925   Float_t mass = pmod/(beta*(Float_t)gamma);
1926
1927   return mass;
1928 }
1929
1930
1931 void AliRICHRecon::FillHistograms()
1932 {
1933
1934   Float_t FittedTrackTheta, FittedTrackPhi;
1935
1936   Float_t ThetaCerenkov    = GetThetaCerenkov();
1937   if(ThetaCerenkov == 999.) return;
1938
1939   Float_t VertZ = GetEventVertexZ();
1940
1941   Float_t TrackTheta = GetTrackTheta();
1942   Float_t TrackPhi   = GetTrackPhi();
1943   Float_t pmod       = GetTrackMomentum();
1944   Float_t pt         = GetTrackPt();
1945   Float_t TrackEta   = GetTrackEta();
1946   Int_t q            = GetTrackCharge();
1947   Float_t TPCLastZ   = GetTrackTPCLastZ();
1948   Float_t MinDist    = GetMinDist();
1949
1950   FittedTrackTheta = GetFittedTrackTheta();
1951   FittedTrackPhi   = GetFittedTrackPhi();
1952   Int_t FittedNPhotonHough = GetFittedHoughPhotons();
1953
1954   if(fDebug)
1955     {
1956       cout << " p " << pmod  << " ThetaC " << ThetaCerenkov
1957            << " rings " << NRings << endl;
1958     }
1959
1960   Int_t NPhotonHough     = GetHoughPhotons();
1961   Float_t NPhotonHoughNorm = GetHoughPhotonsNorm();
1962   Float_t InRing = GetPortionOfRing();
1963
1964   Float_t MassOfParticle = FindMassOfParticle();
1965
1966   Float_t HoughArea = GetHoughArea();
1967   Float_t Multiplicity = GetEventMultiplicity();
1968
1969 //  cout << " area " << HoughArea << " mult " << Multiplicity << endl;
1970
1971   Float_t var[20];
1972
1973 //  var[0] = (Float_t)runID;
1974 //  var[1] = (Float_t)evID;
1975   var[0] = 0;
1976   var[1] = 0;
1977   var[2] = VertZ;
1978   var[3] = pmod;
1979   var[4] = pt;
1980   var[5] = TrackEta;
1981   var[6] = TrackTheta;
1982   var[7] = TrackPhi;
1983   var[8] = FittedTrackTheta;
1984   var[9] = FittedTrackPhi;
1985   var[10] = q;
1986   var[11] = ThetaCerenkov;
1987   var[12] = (Float_t)NPhotonHough;
1988   var[13] = (Float_t)FittedNPhotonHough;
1989   var[14] = InRing;
1990   var[15] = MassOfParticle;
1991   var[16] = HoughArea;
1992   var[17] = Multiplicity;
1993   var[18] = TPCLastZ;
1994   var[19] = MinDist;
1995
1996   hn->Fill(var);
1997
1998   h1_mass->Fill(MassOfParticle);
1999   h2_mvsp->Fill(pmod,MassOfParticle);
2000   h2_mvst->Fill(ThetaCerenkov,MassOfParticle);
2001
2002   FittedTrackTheta = GetFittedTrackTheta();
2003   FittedTrackPhi = GetFittedTrackPhi();
2004
2005   Float_t DiffTheta = FittedTrackTheta - TrackTheta;
2006   Float_t DiffPhi = FittedTrackPhi - TrackPhi;
2007
2008   h1_diffTrackTheta -> Fill(DiffTheta);
2009   h1_diffTrackPhi -> Fill(DiffPhi);
2010
2011   if(ThetaCerenkov > 0.505 && ThetaCerenkov < 0.605)
2012     {
2013       SetPhotonEnergy(6.85);
2014       SetFreonRefractiveIndex();
2015
2016       Float_t pmom = GetTrackMomentum();
2017       Float_t beta = 1./(cos(ThetaCerenkov)*GetFreonRefractiveIndex());
2018       Float_t gamma = 1./sqrt(1.-beta*beta);
2019
2020       Float_t pmomnew = 0.93828*beta*gamma;
2021       Float_t deltap = pmomnew - pmom;
2022       h1_deltap->Fill(deltap);
2023       Float_t deltapop = deltap/pmom;
2024       h1_deltapop->Fill(deltapop);
2025
2026       h1_nprotons->Fill((Float_t)NPhotonHoughNorm);
2027     }
2028
2029   if(q > 0)
2030     {
2031       h2_tvsppos->Fill(pmod,ThetaCerenkov);
2032       hp_1pos->Fill(ThetaCerenkov,(Float_t)NPhotonHough);
2033       hp_1posnorm->Fill(ThetaCerenkov,(Float_t)NPhotonHoughNorm);
2034       h2_1pos->Fill(pmod,(Float_t)NPhotonHough);
2035       h2_1posnorm->Fill(pmod,(Float_t)NPhotonHoughNorm);
2036       h1_houghpos->Fill(ThetaCerenkov);
2037     }
2038 else
2039   {
2040       h2_tvspneg->Fill(pmod,ThetaCerenkov);
2041       hp_1neg->Fill(ThetaCerenkov,(Float_t)NPhotonHough);
2042       hp_1negnorm->Fill(ThetaCerenkov,(Float_t)NPhotonHoughNorm);
2043       h2_1neg->Fill(pmod,(Float_t)NPhotonHough);
2044       h2_1negnorm->Fill(pmod,(Float_t)NPhotonHoughNorm);
2045       h1_houghneg->Fill(ThetaCerenkov);
2046   }
2047
2048   Int_t NPhotons = GetPhotonsNumber();
2049
2050   for (Int_t j=0; j < NPhotons;j++)
2051
2052     {
2053       SetPhotonIndex(j);
2054
2055       Float_t eta = GetPhotonEta();
2056
2057       if(GetPhotonFlag() == 2)
2058         {
2059           h1_photacc->Fill(eta);
2060           Float_t photaccspread = eta - ThetaCerenkov;
2061           h1_photaccspread->Fill(photaccspread);
2062         }
2063
2064     }
2065 }
2066
2067 void AliRICHRecon::Minimization()
2068 {
2069
2070   Double_t arglist;
2071   Int_t ierflag = 0;
2072
2073   static Double_t vstart[2];
2074   static Double_t lower[2], upper[2];
2075   static Double_t step[2]={0.001,0.001};
2076
2077   Double_t TrackThetaNew,TrackPhiNew;
2078   TString chname;
2079   Double_t eps, b1, b2;
2080   Int_t ierflg;
2081
2082   gMyMinuit = new TMinuit(2);
2083   gMyMinuit->SetObjectFit((TObject *)this);
2084   gMyMinuit->SetFCN(fcn);
2085   gMyMinuit->mninit(5,10,7);
2086
2087   vstart[0] = (Double_t)GetTrackTheta();
2088   vstart[1] = (Double_t)GetTrackPhi();
2089
2090   lower[0] = vstart[0] - 0.03;
2091   if(lower[0] < 0) lower[0] = 0.;
2092   upper[0] = vstart[0] + 0.03;
2093   lower[1] = vstart[1] - 0.03;
2094   upper[1] = vstart[1] + 0.03;
2095
2096
2097   gMyMinuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag);
2098   gMyMinuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag);
2099
2100   arglist = -1;
2101
2102   //  gMyMinuit->FixParameter(0);
2103
2104   gMyMinuit->SetPrintLevel(-1);
2105 //  gMyMinuit->mnexcm("SET PRI",&arglist, 1, ierflag);
2106   gMyMinuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
2107   gMyMinuit->mnexcm("SET NOW",&arglist, 1, ierflag);
2108   arglist = 1;
2109   gMyMinuit->mnexcm("SET ERR", &arglist, 1,ierflg);
2110   arglist = -1;
2111
2112   //  gMyMinuit->mnscan();
2113
2114 //  gMyMinuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
2115   gMyMinuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
2116   gMyMinuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
2117
2118   gMyMinuit->mnpout(0,chname, TrackThetaNew, eps , b1, b2, ierflg);
2119   gMyMinuit->mnpout(1,chname, TrackPhiNew, eps , b1, b2, ierflg);
2120
2121   //values after the fit...
2122   SetFittedTrackTheta((Float_t)TrackThetaNew);
2123   SetFittedTrackPhi((Float_t)TrackPhiNew);
2124
2125   delete gMyMinuit;
2126
2127 }
2128
2129 void AliRICHRecon::EstimationOfTheta()
2130 {
2131
2132   Int_t NPhotons = 0;
2133
2134   Float_t ShiftX = GetShiftX();
2135   Float_t ShiftY = GetShiftY();
2136
2137   Float_t *CandidatePhotonX = GetCandidatePhotonX();
2138   Float_t *CandidatePhotonY = GetCandidatePhotonY();
2139
2140   Int_t NPhotonsCandidates = GetCandidatePhotonsNumber();
2141
2142   //  cout << "MINIM: Nphotons " << NPhotonsCandidates << endl;
2143
2144   for (Int_t j=0; j < NPhotonsCandidates; j++)
2145     {
2146
2147       SetPhotonIndex(j);
2148
2149       if(!GetPhotonFlag()) continue;
2150
2151       Float_t Xtoentr = CandidatePhotonX[j] - ShiftX;
2152       Float_t Ytoentr = CandidatePhotonY[j] - ShiftY;
2153
2154       SetEntranceX(Xtoentr);
2155       SetEntranceY(Ytoentr);
2156
2157       FindPhiPoint();
2158
2159       FindThetaPhotonCerenkov();
2160
2161       Float_t ThetaPhotonCerenkov = GetThetaPhotonCerenkov();
2162
2163       //      cout << " ACCEPTED!!! " << ThetaPhotonCerenkov << endl;
2164
2165       SetPhotonEta(ThetaPhotonCerenkov);
2166
2167       NPhotons++;
2168
2169     }
2170
2171   Float_t xmean = 0.;
2172   Float_t x2mean = 0.;
2173   Int_t nev = 0;
2174
2175   for (Int_t j=0; j < NPhotonsCandidates;j++)
2176     {
2177       SetPhotonIndex(j);
2178
2179       Float_t eta = GetPhotonEta();
2180
2181       if(eta != -999.)
2182         {
2183           if(GetPhotonFlag() == 2)
2184             {
2185               xmean += eta;
2186               x2mean += eta*eta;
2187               nev++;
2188             }
2189         }
2190     }
2191
2192   if(nev > 0)
2193     {
2194       xmean /=(Float_t)nev;
2195       x2mean /=(Float_t)nev;
2196     } else {
2197       xmean = 0.;
2198       x2mean = 0.;
2199     }
2200
2201   Float_t RMS = sqrt(x2mean - xmean*xmean);
2202
2203   //  cout << " RMS " << RMS;
2204
2205   SetEstimationOfTheta(xmean);
2206   SetEstimationOfThetaRMS(RMS);
2207 }
2208
2209 void fcn(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t iflag)
2210 {
2211   AliRICHRecon *gMyRecon = (AliRICHRecon*)gMyMinuit->GetObjectFit();
2212
2213   Float_t p0 = (Float_t)par[0];
2214   Float_t p1 = (Float_t)par[1];
2215
2216   gMyRecon->SetTrackTheta(p0);
2217   gMyRecon->SetTrackPhi(p1);
2218
2219   gMyRecon->EstimationOfTheta();
2220   Float_t RMS = gMyRecon->GetEstimationOfThetaRMS();
2221
2222   Int_t HoughPhotons = gMyRecon->GetHoughPhotons();
2223
2224
2225   f = (Double_t)(1000*RMS/(Float_t)HoughPhotons);
2226
2227   if(fDebug) cout << "   f   " << f
2228                   << " theta " << par[0] << " phi " << par[1]
2229                   << " HoughPhotons " << HoughPhotons << endl;
2230
2231   if(fDebug&&iflag == 3)
2232     {
2233             cout << " --- end convergence...summary --- " << endl;
2234             cout << " theta " << par[0] << endl;
2235             cout << "  phi  " << par[1] << endl;
2236     }
2237 }
2238
2239 void AliRICHRecon::waiting()
2240 {
2241   if(!kDISPLAY) return;
2242   cout << " Press any key to continue...";
2243
2244 //  gSystem->ProcessEvents();
2245   getchar();
2246
2247   cout << endl;
2248
2249   return;
2250 }
2251
2252 /*
2253 void ~AliRICHRecon()
2254 {
2255 }
2256 */