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