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