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