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