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