]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/AliRICH.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Do not export any Make-depend, so we had to change the name of the
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
CommitLineData
fe4da5cc 1///////////////////////////////////////////////////////////////////////////////
2// //
3// Ring Imaging Cherenkov //
4// This class contains the basic functions for the Ring Imaging Cherenkov //
5// detector. Functions specific to one particular geometry are //
6// contained in the derived classes //
7// //
8//Begin_Html
9/*
10<img src="gif/AliRICHClass.gif">
11*/
12//End_Html
13// //
14// //
15///////////////////////////////////////////////////////////////////////////////
16
17#include <TBRIK.h>
18#include <TNode.h>
19#include <TRandom.h>
20
21#include <TClass.h>
22#include "AliRICH.h"
23#include "AliRun.h"
24#include "TGeant3.h"
25
26
27ClassImp(AliRICH)
28
29//_____________________________________________________________________________
30AliRICH::AliRICH()
31{
32 //
33 // Default constructor for RICH
34 //
35 fIshunt = 0;
36 fHits = 0;
37 fMips = 0;
38 fCkovs = 0;
39 fPadhits = 0;
40 fNmips = 0;
41 fNckovs = 0;
42 fNpadhits = 0;
43
44 fChslope = 0;
45 fAlphaFeed= 0;
46 fSxcharge = 0;
47 fIritri = 0;
48}
49
50//_____________________________________________________________________________
51AliRICH::AliRICH(const char *name, const char *title)
52 : AliDetector(name,title)
53{
54 //
55 // Standard constructor for RICH
56 //
57 fIshunt = 0;
58 fNmips = 0;
59 fNckovs = 0;
60 fNpadhits = 0;
61 //
62 // Allocate space for different components
63 fHits = new TClonesArray("AliRICHhit", 100);
64 fMips = new TClonesArray("AliRICHmip", 100);
65 fCkovs = new TClonesArray("AliRICHckov", 100);
66 fPadhits = new TClonesArray("AliRICHpadhit", 100);
67 //
68 // Set parameters to default value
69 fChslope = 40;
70 fAlphaFeed= 0.04;
71 fSxcharge = 0.18;
72 fIritri = 0;
73 //
74 SetMarkerColor(6);
75 SetMarkerStyle(20);
76 SetMarkerSize(0.5);
77}
78
79//_____________________________________________________________________________
80AliRICH::~AliRICH()
81{
82 //
83 // Destructor for RICH
84 //
85 fIshunt = 0;
86 delete fHits;
87 fMips->Delete(); delete fMips;
88 fCkovs->Delete(); delete fCkovs;
89 fPadhits->Delete(); delete fPadhits;
90}
91
92//_____________________________________________________________________________
93void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
94{
95 //
96 // Add a RICH hit
97 //
98 switch ((int) hits[0]) {
99 case 0:
100 {
101 //
102 // Simple hit
103 TClonesArray &lhits = *fHits;
104 new(lhits[fNhits++]) AliHit(fIshunt,track);
105 AliHit *lhit = (AliHit*) fHits->AddrAt(fNhits-1);
106 lhit->fX = hits[19];
107 lhit->fY = hits[20];
108 lhit->fZ = hits[21];
109 } break;
110 case 1:
111 //
112 // MIP hit
113 AddMipHit(track,vol,hits);
114 break;
115 case 2:
116 //
117 // Cherenkov hit
118 AddCkovHit(track,vol,hits);
119 break;
120 case 3:
121 //
122 // Pad hit
123 AddPadHit(track,vol,hits);
124 break;
125 case 4:
126 // Update a mip hit
127 UpdateMipHit(hits);
128 break;
129 default:
130 printf("Error: AliRICH::AddHit flag %d not defined./n",(int) hits[0]);
131 return;
132 }
133}
134
135//_____________________________________________________________________________
136void AliRICH::AddMipHit(Int_t track, Int_t *vol, Float_t *hits)
137{
138 // Adds a mip hit in the RICH.
139 TClonesArray &lhits = *fMips;
140 new(lhits[fNmips++]) AliRICHmip(fIshunt,track,vol,hits,
141 fNckovs,fNpadhits);
142}
143
144//_____________________________________________________________________________
145void AliRICH::AddCkovHit(Int_t track, Int_t *vol, Float_t *hits)
146{
147 //
148 // Adds a cerenkov hit in the RICH.
149 //
150 TClonesArray &lhits = *fCkovs;
151 AliRICHmip *lmip = (AliRICHmip*) fMips->AddrAt(fNmips-1);
152 //
153 // If this ckov come from a mip update the mip lastckov entry.
154 Int_t fmipslocal=-1;
155 if (lmip->GetZ() != -999.0) {
156 fmipslocal = fNmips-1;
157 lmip->SetLastCkov(fNckovs);
158 }
159 new(lhits[fNckovs++]) AliRICHckov(fIshunt,track,vol,hits,
160 fmipslocal,fNpadhits);
161}
162
163//_____________________________________________________________________________
164void AliRICH::AddPadHit(Int_t track, Int_t *vol, Float_t *hits)
165{
166 //
167 // Adds pad hits in the RICH
168 //
169 TClonesArray &lhits = *fPadhits;
170 // Update the last padhit of the respective particle:
171 if ((int) hits[1]==50) { // a ckov
172 ((AliRICHckov *) fCkovs->AddrAt(fNckovs-1))->SetLastpad(fNpadhits);
173 new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,-1,fNckovs-1);
174 }else { // a mip
175 ((AliRICHmip *) fMips->AddrAt(fNmips-1))->SetLastpad(fNpadhits);
176 new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,fNmips-1,-1);
177 }
178}
179
180//_____________________________________________________________________________
181void AliRICH::BuildGeometry()
182{
183 //
184 // Builds a TNode geometry for event display
185 //
186 TNode *Node, *Top;
187
188 const int kColorRICH = kGreen;
189 //
190 Top=gAlice->GetGeometry()->GetNode("alice");
191
192 new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
193 new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
194 new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
195 new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
196 new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
197 new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
198 new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
199 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
200 Top->cd();
201 Node = new TNode("RICH1","RICH1","S_RICH",0,471.8999,165.2599,"rot993");
202 Node->SetLineColor(kColorRICH);
203 fNodes->Add(Node);
204 Top->cd();
205 Node = new TNode("RICH2","RICH2","S_RICH",171,470,0,"rot994");
206 Node->SetLineColor(kColorRICH);
207 fNodes->Add(Node);
208 Top->cd();
209 Node = new TNode("RICH3","RICH3","S_RICH",0,500,0,"rot995");
210 Node->SetLineColor(kColorRICH);
211 fNodes->Add(Node);
212 Top->cd();
213 Node = new TNode("RICH4","RICH4","S_RICH",-171,470,0,"rot996");
214 Node->SetLineColor(kColorRICH);
215 fNodes->Add(Node);
216 Top->cd();
217 Node = new TNode("RICH5","RICH5","S_RICH",161.3999,443.3999,-165.3,"rot997");
218 Node->SetLineColor(kColorRICH);
219 fNodes->Add(Node);
220 Top->cd();
221 Node = new TNode("RICH6","RICH6","S_RICH",0,471.8999,-165.3,"rot998");
222 Node->SetLineColor(kColorRICH);
223 fNodes->Add(Node);
224 Top->cd();
225 Node = new TNode("RICH7","RICH7","S_RICH",-161.399,443.3999,-165.3,"rot999");
226 Node->SetLineColor(kColorRICH);
227 fNodes->Add(Node);
228}
229
230//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
231//
232// Section for Rich Step
233#define MAXPH 1000
234
235static Float_t sVloc[3];
236static Float_t sVectIn[3];
237static Float_t sDetot;
238static Float_t sYanode[MAXPH];
239static Float_t sXpad[MAXPH];
240static Float_t sYpad[MAXPH];
241static Int_t sNpx;
242static Int_t sNpy;
243static Float_t sDxp;
244static Float_t sDyp;
245static Float_t sDlx;
246static Float_t sDly;
247static Float_t sDpad;
248
249
250static Int_t sNsecon;
251//static Float_t sQint[2];
252
253#define MAXSEC 2000
254
255static Int_t sNpads;
256static Int_t sIpx[MAXSEC];
257static Int_t sIpy[MAXSEC];
258static Int_t sIqpad[MAXSEC];
259static Int_t sNphlink[MAXSEC];
260static Int_t sNphoton;
261
262static Int_t sNfeed, sNfeedd, sNdir;
263
264static Float_t sPparent;
265static Float_t sThincParent;
266static Int_t sIloss[MAXPH];
267static Int_t sIprod[MAXPH];
268static Float_t sXphit[MAXPH];
269static Float_t sYphit[MAXPH];
270
271static Float_t sSycharge;
272static Float_t sXsox, sYsox, sZsox;
273
274static Int_t sMckov[MAXPH];
275static Int_t idpartgx;
276static Float_t phisx;
277static Int_t nmodsx;
278static Float_t psx;
279static Float_t xsx;
280static Float_t ysx;
281static Float_t thetasx;
282
283
284//_____________________________________________________________________________
285void AliRICH::Init()
286{
287 //
288 // Initialise RICH after that it has been built
289 //
290 const Float_t sconv=2*TMath::Sqrt(2*TMath::Log(2));
291 const Float_t yK3=1.20;
292 Float_t ansp;
293 Int_t i;
294 //
295 sNpx=162;
296 sNpy=162;
297 sDxp=0.80;
298 sDyp=0.80;
299 ansp=sDyp/2;
300 sDlx=sNpx*sDxp/2;
301 sDly=sNpy*sDyp/2;
302 sDpad=0.2;
303 //
304 for(i=0;i<sNpx;i++) {
305 sXpad[i]=i*sDxp;
306 sYpad[i]=i*sDyp;
307 }
308 for(i=0;i<2*sNpy+1;i++) sYanode[i]=ansp/2+i*ansp;
309 //
310
311 sSycharge=4*TMath::ATanH(1/TMath::Sqrt(2*yK3))/TMath::Pi()
312 /(1-0.5*TMath::Sqrt(yK3))/sDpad/sconv;
313 //
314}
315
316
317//___________________________________________________________________________
318void AliRICH::StepManager()
319{
320 //
321 // Called at every step in the RICH
322 //
323
324 AliMC* pMC = AliMC::GetMC();
325 TGeant3 *geant3 = (TGeant3*) pMC;
326
327 const Float_t xshift[3] = { 41.3, 0, -41.3 };
328 static Float_t polar[3] = {0, 0, 0};
329 const Int_t nrooth = 25;
330
331 static Int_t ixold=-1, iyold=-1;
332
333 // System generated locals
334 Int_t j, i1;
335 Float_t r1, r2;
336
337 // Local variables
338 Float_t ranf[2], rrhh[nrooth], phiangle, cost, vmod;
339 //Int_t idpartsx;
340 Int_t i;
341 Float_t t, vxloc[3];
342 Int_t ll, irivol[2];
343 Int_t lcase;
344 //Int_t iprimx;
345 Int_t ix, iy;
346 Float_t stwght;
347 Int_t ncher;
348 Float_t cophi;
349 Float_t dir[3];
350 Int_t ihitrak;
351 Int_t medprod;
352
353 Int_t nmult=0;
354 //Float_t xtrig[200], ytrig[200];
355 //Int_t itrig[200];
356
357
358
359 // ILOSS = 0 NOT LOST
360 // 1 REFLECTED FREON-QUARZ
361 // 2 REFLECTED QUARZ-METHANE
362 // 3 REFLECTED METHANE-CSI
363 // 11 ABSORBED IN FREON
364 // 12 ABSORBED IN QUARZ
365 // 13 ABSORBED IN METHANE
366 // 21 CSI QE
367 // IPROD = 1 PRODUCED IN FREON
368 // IPROD = 2 PRODUCED IN QUARZ
369
370 // new (changed NROOTH from 10 to 25!!!!!!!!!!!!!)
371
372
373 Int_t *idtmed = gAlice->Idtmed();
374
375 //--------------------------------------------------------------------------
376
377 // MIP inside CsI
378
379 if (geant3->Gckine()->charge) {
380
381 // Charged particles treatment
382 if (fIritri && !geant3->Gctrak()->upwght) {
383 if (geant3->Gctmed()->numed == idtmed[fIritri-1]) {
384 if (geant3->Gcking()->ngkine > 0) {
385 strncpy((char *)&lcase,"HADR",4);
386 if (geant3->Gcking()->kcase == lcase) {
387 i1 = geant3->Gcking()->ngkine;
388 for (i = 1; i <= i1; ++i) {
389 pMC->Gmtod(geant3->Gckin3()->gpos[i-1], vxloc, 1);
390 pMC->Gmtod(geant3->Gcking()->gkin[i-1], dir, 2);
391 if (geant3->Gcking()->gkin[i-1][4] == 8. ||
392 geant3->Gcking()->gkin[i-1][4] == 9.) {
393 ++nmult;
394 // Computing 2nd power
395 r1 = dir[0];
396 // Computing 2nd power
397 r2 = dir[2];
398 //theta = TMath::ATan2(TMath::Sqrt(r1*r1+r2*r2),dir[1]);
399 //xtrig[nmult - 1] = theta;
400 //ytrig[nmult - 1] = vxloc[1] + .25;
401 //itrig[nmult - 1] = (Int_t) geant3->Gcking()->gkin[i-1][4];
402 }
403 }
404 }
405 }
406 }
407 if ((geant3->Gctmed()->numed == idtmed[1006-1] &&
408 geant3->Gctrak()->inwvol == 2) ||
409 geant3->Gctrak()->istop) {
410 if (!nmult) {
411 printf("NOT TRIGGERED\n");
412 sDetot = 0.;
413 sNsecon = 0;
414 sNpads = 0;
415 sNphoton = 0;
416 sNfeed = 0;
417 sNfeedd = 0;
418 sNdir = 0;
419 geant3->Gctrak()->istory = 0;
420 geant3->Gctrak()->upwght = 0.;
421 geant3->Gcflag()->ieotri = 1;
422 nmult = 0;
423 //sQint[0] = 0.;
424 //sQint[1] = 0.;
425 } else {
426 printf("TRIGGERED %d\n",nmult);
427 }
428 }
429 }
430 // MIP inside Methane
431 if (geant3->Gctmed()->numed == idtmed[1009-1]) {
432
433 // new
434 // If particle produced already Cerenkov Photons (istory=1)
435 // update the impact point only
436 if (geant3->Gctrak()->istory == 1) {
437 // Direction of incidence and where did it hit ?
438 pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
439 pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
440 phiangle = TMath::ATan2(dir[2], dir[0]);
441 if (phiangle < 0.) phiangle += 2*TMath::Pi();
442 i1 = nrooth;
443 for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
444 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
445 irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
446 // NMODSX
447 ihitrak = gAlice->CurrentTrack();
448 rrhh[0] = 4.;
449 // flag to say this is update
450 rrhh[1] = sVloc[0] + sDlx;
451 // XSX
452 rrhh[2] = sVloc[2] + sDly;
453 // YSX
454 rrhh[3] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
455 // NMODSX
456 // Computing 2nd power
457 r1 = dir[0];
458 // Computing 2nd power
459 r2 = dir[2];
460 rrhh[4] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
461 // theta
462 rrhh[5] = phiangle;
463 // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
464 AddHit(ihitrak,irivol,rrhh);
465 }
466 // enew
467 // Record particle properties
468 // If particle produced already Cerenkov Photons (istory=1)
469
470 // update the impact point only
471 if (geant3->Gctrak()->istory != 2) {
472 if (!geant3->Gctrak()->istory) {
473 ++sNsecon;
474
475 // Is this a primary particle ?
476 //iprimx = 1;
477 //if (geant3->Gctrak()->upwght) iprimx = 0;
478
479 // Where did it come from ?
480 sXsox = geant3->Gckine()->vert[0];
481 sYsox = geant3->Gckine()->vert[1];
482 sZsox = geant3->Gckine()->vert[2];
483
484 // Momentum
485 psx = geant3->Gctrak()->vect[6];
486
487 // Particle type and parent
488 //idpartsx = geant3->Gckine()->ipart;
489 r1 = geant3->Gctrak()->upwght / 100.;
490 idpartgx = Int_t(r1+0.5);
491 if (!geant3->Gctrak()->upwght) {
492 sPparent = geant3->Gctrak()->vect[6];
493 sThincParent = thetasx;
494 }
495
496 // Direction of incidence and where did it hit ?
497 pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
498 pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
499 // Computing 2nd power
500 r1 = dir[0];
501 // Computing 2nd power
502 r2 = dir[2];
503 thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
504 phisx = TMath::ATan2(dir[2], dir[0]);
505 if (phisx < 0.) phisx += 2*TMath::Pi();
506 ysx = sVloc[2] + sDly;
507 xsx = sVloc[0] + sDlx;
508 nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
509 // new
510 i1 = nrooth;
511 for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
512 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
513 irivol[1] = nmodsx;
514 ihitrak = gAlice->CurrentTrack();
515 rrhh[0] = 1.;
516 // Flag to say this is MIP
517 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
518 rrhh[2] = xsx;
519 rrhh[3] = ysx;
520 rrhh[4] = (Float_t) nmodsx;
521 // Module Number
522 rrhh[5] = thetasx;
523 rrhh[6] = geant3->Gctrak()->tofg;
524 // in seconds
525 rrhh[7] = (Float_t) idpartgx;
526 // mips specific
527 rrhh[8] = phisx;
528 rrhh[9] = psx;
529 // charge of current particle in electron charge unit;
530 rrhh[10] = geant3->Gckine()->charge;
531 rrhh[11] = -999.;
532 // Zo of ckov generation
533 rrhh[12] = 0.;
534 // no ckov !!!
535 AddHit(ihitrak, irivol, rrhh);
536 // end of new
537
538 // Earmark track as being recorded in methane gap
539
540 geant3->Gctrak()->istory = 2;
541 }
542 }
543
544 // Signal generation in methane gap
545 pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
546 pMC->Gmtod(geant3->Gckine()->vert, vxloc, 1);
547 ix = (Int_t) ((sVloc[0] + sDlx) / sDxp);
548 iy = (Int_t) ((sVloc[2] + sDly) / sDyp);
549
550 // Is this the first step?
551 if (ixold == -1 && iyold == -1) {
552 ixold = ix;
553 iyold = iy;
554 for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
555 }
556
557 // Mip left gap
558 if (geant3->Gctrak()->inwvol == 2 || geant3->Gctrak()->istop) {
559 sDetot += geant3->Gctrak()->destep;
560 if (sDetot > 0.) RichIntegration();
561 sDetot = 0.;
562 ixold = -1;
563 iyold = -1;
564
565 // Mip left current pad
566 } else if (ixold != ix || iyold != iy) {
567 if (sDetot > 0.) RichIntegration();
568 for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
569 sDetot = geant3->Gctrak()->destep;
570 ixold = ix;
571 iyold = iy;
572 } else {
573 sDetot += geant3->Gctrak()->destep;
574 }
575 }
576 }
577
578 // End charged particles treatment
579 //
580 // Treat photons produced in Freon and Quartz
581 if (geant3->Gckin2()->ngphot > 0 &&
582 (geant3->Gctmed()->numed == idtmed[1004-1] ||
583 geant3->Gctmed()->numed == idtmed[1003-1])) {
584 if (!geant3->Gctrak()->upwght) {
585
586 // If it is a primary, save all generated photons
587 i1 = geant3->Gckin2()->ngphot;
588 for (i = 1; i <= i1; ++i) {
589 ++sNphoton;
590 if (sNphoton > MAXPH) {
591 sNphoton = MAXPH;
592 printf("ATTENTION NPHOTON %d\n",sNphoton);
593 continue;
594 }
595
596 // Production medium
597 medprod = 1;
598 if (geant3->Gctmed()->numed == idtmed[1003-1]) medprod = 2;
599 //
600 // Production angle and energy
601 vmod=0;
602 cost=0;
603 for(j=0;j<3;j++) {
604 cost+=geant3->Gckin2()->xphot[i-1][3+j]*geant3->Gctrak()->vect[3+j];
605 vmod+=geant3->Gckin2()->xphot[i-1][3+j]*
606 geant3->Gckin2()->xphot[i-1][3+j];
607 }
608 cost/=sqrt(vmod);
609 sIloss[sNphoton - 1] = 22;
610 sIprod[sNphoton - 1] = medprod;
611 sXphit[sNphoton - 1] = 0.;
612 sYphit[sNphoton - 1] = 0.;
613 stwght = geant3->Gctrak()->upwght;
614 geant3->Gctrak()->upwght = (Float_t) sNphoton;
615 geant3->Gskpho(i);
616 gAlice->SetTrack(0, gAlice->CurrentTrack(), 50,
617 &geant3->Gckin2()->xphot[i-1][3],geant3->Gckin2()->xphot[i-1],
618 polar,geant3->Gctrak()->tofg,"Cherenkov", ncher);
619 sMckov[sNphoton - 1] = ncher;
620 geant3->Gctrak()->upwght = stwght;
621 }
622 } else {
623 stwght = geant3->Gctrak()->upwght;
624 geant3->Gctrak()->upwght = 0.;
625 geant3->Gskpho(0);
626 geant3->Gctrak()->upwght = stwght;
627 }
628
629 // Particle did not yet pass the methane gap
630 if (geant3->Gctrak()->istory == 0) {
631 geant3->Gctrak()->istory = 1;
632 ++sNsecon;
633 // Is this a primary particle ?
634 //iprimx = 1;
635 //if (geant3->Gctrak()->upwght) iprimx = 0;
636
637 // Where did it come from ?
638 sXsox = geant3->Gckine()->vert[0];
639 sYsox = geant3->Gckine()->vert[1];
640 sZsox = geant3->Gckine()->vert[2];
641
642 // Where did it hit ?
643 pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
644 pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
645 ysx = sVloc[2] + sDly;
646 if (geant3->Gctmed()->numed == idtmed[1004-1]) {
647 nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
648 xsx = sVloc[0] + sDlx +
649 xshift[geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 2] - 1];
650 } else if (geant3->Gctmed()->numed == idtmed[1003-1]) {
651 nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
652 xsx = sVloc[0] + sDlx;
653 } else {
654 nmodsx = 0;
655 }
656
657 // Momentum and direction of incidence
658 psx = geant3->Gctrak()->vect[6];
659 // Computing 2nd power
660 r1 = dir[0];
661 // Computing 2nd power
662 r2 = dir[2];
663 thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
664 phisx = TMath::ATan2(dir[2], dir[0]);
665 if (phisx < 0.) phisx += 2*TMath::Pi();
666
667 // Particle type and parent
668 //idpartsx = geant3->Gckine()->ipart;
669 r1 = geant3->Gctrak()->upwght / 100.;
670 idpartgx = Int_t(r1+0.5);
671 if (!geant3->Gctrak()->upwght) {
672 sPparent = geant3->Gctrak()->vect[6];
673 sThincParent = thetasx;
674 }
675 // new
676 for (ll = 0; ll <nrooth; ++ll) rrhh[ll] = 0;
677 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
678 irivol[1] = nmodsx;
679 ihitrak = gAlice->CurrentTrack();
680 rrhh[0] = 1.;
681 // Flag to say that this is MIP
682 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
683 rrhh[2] = xsx;
684 rrhh[3] = ysx;
685 rrhh[4] = (Float_t) nmodsx;
686 // Module Number
687 rrhh[5] = thetasx;
688 rrhh[6] = geant3->Gctrak()->tofg;
689 // in seconds
690 rrhh[7] = (Float_t) idpartgx;
691 // mips specific
692 rrhh[8] = phisx;
693 rrhh[9] = psx;
694 // charge of current particle in electron charge unit;
695 rrhh[10] = geant3->Gckine()->charge;
696 rrhh[11] = sVloc[1];
697 // Zo of ckov generation
698 rrhh[12] = 1.;
699 // ckov generation
700 AddHit(ihitrak, irivol, rrhh);
701 // enew
702 }
703 }
704
705 // Current particle is cherenkov photon
706 if (geant3->Gckine()->ipart == 50) {
707 pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
708 // WRITE(6,* ) UPWGHT, VLOC(2), NUMED, DESTEP
709 // Photon crosses ch4-csi boundary
710 // take into account fresnel losses with complex refraction index
711 if (geant3->Gctrak()->inwvol == 1 && geant3->Gctmed()->numed == idtmed[1006-1]) {
712
713 // fresnel losses commented out for the moment
714 // make sure that qe correction for fresnel losses is also switched off !
715 // CALL FRESNELCSI
716 // IF (ISTOP .EQ. 2) RETURN
717 // Put transmission of electrodes in by hand
718 pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
719 cophi = TMath::Cos(TMath::ATan2(dir[0], dir[1]));
720 t = (1. - .025 / cophi) * (1. - .05 / cophi);
721 pMC->Rndm(ranf, 1);
722 if (ranf[0] > t) {
723 if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5)<MAXPH)
724 sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 15;
725 geant3->Gctrak()->istop = 2;
726 return;
727 }
728 }
729
730 // Photon lost energy in CsI
731 if (geant3->Gctrak()->destep > 0. && geant3->Gctmed()->numed == idtmed[1006-1]) {
732 geant3->Gctrak()->istop = 2;
733 r1 = geant3->Gctrak()->upwght / 100.;
734 if (Int_t(r1+0.5) > 50) {
735 ++sNfeedd;
736 } else {
737 ++sNdir;
738 }
739 // WRITE(6,*) 'PHOTON',UPWGHT, MAXPH
740 if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH)
741 sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 0;
742 // write(6,*) 'photon detected'
743 for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
744 // new
745 // copied from miphit in Freon or Quartz
746 // Where did it hit ?
747 pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
748
749 // Momentum and direction of incidence
750 for (ll = 0; ll < nrooth; ++ll) rrhh[ll]=0;
751 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
752 // ???
753 irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
754 // ???
755 ihitrak = gAlice->CurrentTrack();
756 rrhh[0] = 2.;
757 // Flag to say that this is CK
758 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
759 rrhh[2] = sVloc[0] + sDlx;
760 rrhh[3] = sVloc[2] + sDly;
761 rrhh[4] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
762 // ??? Module Number
763 // Computing 2nd power
764 r1 = dir[0];
765 // Computing 2nd power
766 r2 = dir[2];
767 rrhh[5] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
768 // THETASX
769 rrhh[6] = geant3->Gctrak()->tofg;
770 // in seconds
771 r1 = geant3->Gctrak()->upwght / 100.;
772 rrhh[7] = (Float_t) Int_t(r1+0.5);
773 // ckov specific
774 // Feedback ???
775 rrhh[8] = geant3->Gctrak()->getot;
776 rrhh[9] = 0.;
777 // Stop in CsI
778 AddHit(ihitrak, irivol, rrhh);
779 // end of new
780 RichIntegration();
781 return;
782 }
783 if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH)
784 {
785 // Losses by reflection and absorption and leaving detector
786 if (sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] == 22) {
787 i1 = geant3->Gctrak()->nmec;
788 for (i = 0; i < i1; ++i) {
789 // Reflection loss
790 if (geant3->Gctrak()->lmec[i] == 106) {
791 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=10;
792 if (geant3->Gctmed()->numed == idtmed[1004-1])
793 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=1;
794
795 if (geant3->Gctmed()->numed == idtmed[1003-1])
796 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=2;
797
798 // Absorption loss
799 } else if (geant3->Gctrak()->lmec[i] == 101) {
800 geant3->Gctrak()->istop = 2;
801 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=20;
802 if (geant3->Gctmed()->numed == idtmed[1004-1])
803 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=11;
804
805 if (geant3->Gctmed()->numed == idtmed[1003-1])
806 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=12;
807
808 if (geant3->Gctmed()->numed == idtmed[1005-1])
809 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
810
811 if (geant3->Gctmed()->numed == idtmed[1009-1])
812 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
813
814 if (geant3->Gctmed()->numed == idtmed[1001-1])
815 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=14;
816
817 // CsI inefficiency
818 if (geant3->Gctmed()->numed == idtmed[1006-1])
819 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=16;
820
821
822 // Photon goes out of tracking scope
823 } else if (geant3->Gctrak()->lmec[i] == 30)
824 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=21;
825
826 }
827 }
828 }
829 }
830}
831
832
833//_____________________________________________________________________________
834void AliRICH::RichIntegration()
835{
836 //
837 // Integrates over RICH pads
838 //
839
840 AliMC* pMC = AliMC::GetMC();
841 TGeant3 *geant3 = (TGeant3*) pMC;
842
843 Int_t i1, i2;
844 Float_t r1;
845
846 // Local variables
847 Float_t rrhh[25];
848 Float_t qtot;
849 Int_t ifeed;
850 Float_t x0;
851 Int_t ixmin, ixmax, iymin, iymax;
852 Int_t ll, ix, iy;
853 Float_t qp;
854 Float_t source[3];
855 Int_t irivol[2];
856 Float_t y0a, qtot_check, xi1, xi2, yi1, yi2;
857 Int_t nph = 0, iqp;
858 Int_t ihitrak, modulen;
859 //Int_t isec[MAXSEC];
860 // ILOSS = 0 NOT LOST
861 // 1 REFLECTED FREON-QUARZ
862 // 2 REFLECTED QUARZ-METHANE
863 // 3 REFLECTED METHANE-CSI
864 // 11 ABSORBED IN FREON
865 // 12 ABSORBED IN QUARZ
866 // 13 ABSORBED IN METHANE
867 // 21 CSI QE
868 // IPROD = 1 PRODUCED IN FREON
869 // IPROD = 2 PRODUCED IN QUARZ
870 // get charge for MIP or cherenkov photon
871
872 if (geant3->Gckine()->ipart == 50) {
873 GetCharge(qtot);
874 modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
875 //sQint[1] = qtot;
876 } else {
877 GetChargeMip(qtot);
878 modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
879 //sQint[0] = qtot;
880 }
881 //
882 pMC->Gmtod(sVectIn, sVloc, 1);
883 if (TMath::Abs(sVloc[0]) >= sDlx) return;
884
885 if (TMath::Abs(sVloc[2]) >= sDly) return;
886
887 sVloc[0] += sDlx;
888 x0 = sVloc[0];
889 sVloc[2] += sDly;
890 //y0 = sVloc[2];
891 AnodicWires(y0a);
892 ixmin = (Int_t) ((x0 - fSxcharge * 5.) / sDxp) + 1;
893 if (ixmin < 1) ixmin = 1;
894 ixmax = (Int_t) ((x0 + fSxcharge * 5.) / sDxp) + 1;
895 if (ixmax > sNpx) ixmax = sNpx;
896 iymin = (Int_t) ((y0a - sSycharge * 5.) / sDyp) + 1;
897 if (iymin < 1) iymin = 1;
898 iymax = (Int_t) ((y0a + sSycharge * 5.) / sDyp) + 1;
899 if (iymax > sNpy) iymax = sNpy;
900 qtot_check = 0.;
901 i1 = ixmax;
902 for (ix = ixmin; ix <= i1; ++ix) {
903 i2 = iymax;
904 for (iy = iymin; iy <= i2; ++iy) {
905 xi1 = (sXpad[ix - 1] - x0) / sDpad;
906 xi2 = xi1 + sDxp / sDpad;
907 yi1 = (sYpad[iy - 1] - y0a) / sDpad;
908 yi2 = yi1 + sDyp / sDpad;
909 qp = qtot * FMathieson(xi1, xi2) * FMathieson(yi1, yi2);
910 iqp = (Int_t) TMath::Abs(qp);
911 qtot_check += TMath::Abs(qp);
912
913 // FILL COMMON FOR PADS
914
915 if (iqp) {
916 if (iqp > 4095) {
917 iqp = 4096;
918 }
919 ++sNpads;
920 if (sNpads > MAXSEC) {
921 // write(6,*) 'attention npads',npads
922 sNpads = MAXSEC;
923 }
924 sIpx[sNpads - 1] = ix;
925 sIpy[sNpads - 1] = iy;
926 sIqpad[sNpads - 1] = iqp;
927 // TAG THE ORIGIN OF THE CHARGE DEPOSITION
928 r1 = geant3->Gctrak()->upwght / 100.;
929 ifeed = Int_t(r1+0.5);
930 sNphlink[sNpads - 1] = 0;
931 if (geant3->Gckine()->ipart != 50) {
932 // MIP
933 //isec[sNpads - 1] = sNsecon;
934 } else {
935 if (ifeed < 50) {
936 // CERENKOV PHOTON
937
938 nph = Int_t(geant3->Gctrak()->upwght+0.5);
939 sNphlink[sNpads - 1] = nph;
940 sXphit[nph - 1] = sVloc[0];
941 sYphit[nph - 1] = sVloc[2];
942 //isec[sNpads - 1] = sNsecon + 100;
943 } else if (ifeed == 51) {
944 // FEEDBACK FROM CERENKOV
945
946 //isec[sNpads - 1] = sNsecon + 300;
947 } else if (ifeed == 52) {
948 // FEEDBACK FROM MIP
949
950 //isec[sNpads - 1] = sNsecon + 200;
951 }
952 }
953 // Generate the hit for Root IO
954 for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
955 irivol[0] = modulen;
956 irivol[1] = nmodsx;
957 rrhh[0] = 0.;
958 // Flag to say this is a hit
959 rrhh[1] = xsx;
960 rrhh[2] = ysx;
961 rrhh[3] = psx;
962 rrhh[4] = thetasx;
963 rrhh[5] = phisx;
964 rrhh[6] = (Float_t) idpartgx;
965 rrhh[7] = sXsox;
966 rrhh[8] = sYsox;
967 rrhh[9] = sZsox;
968 rrhh[10] = (Float_t) sIpx[sNpads - 1];
969 rrhh[11] = (Float_t) sIpy[sNpads - 1];
970 rrhh[12] = (Float_t) sIqpad[sNpads - 1];
971 ihitrak = gAlice->CurrentTrack();
972 if (sNphlink[sNpads - 1] > 0) {
973 rrhh[13] = sPparent;
974 rrhh[14] = sThincParent;
975 rrhh[15] = (Float_t) sIloss[nph - 1];
976 rrhh[16] = (Float_t) sIprod[nph - 1];
977 rrhh[17] = sXphit[nph - 1];
978 rrhh[18] = sYphit[nph - 1];
979 ihitrak = sMckov[nph - 1];
980 }
981 // This is the current position of the hit
982 rrhh[19] = geant3->Gctrak()->vect[0];
983 rrhh[20] = geant3->Gctrak()->vect[1];
984 rrhh[21] = geant3->Gctrak()->vect[2];
985 // PRINT *,' '
986 // PRINT *,'RXAHIT(oldhit)=[',RRHH(20),',',RRHH(21),',',
987 // + RRHH(22),']'
988 // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
989 AddHit(ihitrak, irivol, rrhh);
990 // new Padhits
991 for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
992 rrhh[0] = 3.;
993 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
994 rrhh[2] = (Float_t) sIpx[sNpads - 1];
995 rrhh[3] = (Float_t) sIpy[sNpads - 1];
996 rrhh[4] = (Float_t) modulen;
997 rrhh[5] = -1.;
998 // !!!Prod
999 rrhh[6] = (Float_t) sIqpad[sNpads - 1];
1000 AddHit(ihitrak, irivol, rrhh);
1001 // enew
1002 }
1003 }
1004 }
1005 // IF (IPART .NE. 50) WRITE(6,*) 'CHECK',QTOT,QTOT_CHECK
1006
1007 // GENERATE FEEDBACK PHOTONS
1008
1009 source[0] = x0 - sDlx;
1010 source[1] = sVloc[1] - .2;
1011 source[2] = y0a - sDly;
1012 pMC->Gdtom(source, source, 1);
1013 FeedBack(source, qtot);
1014 return;
1015}
1016
1017//_____________________________________________________________________________
1018void AliRICH::AnodicWires(Float_t &y0a)
1019{
1020 //
1021 // Position of anodic wires
1022 //
1023 Int_t i1;
1024
1025 // Local variables
1026 Int_t i;
1027 Float_t ass_i, y0, ass_i1;
1028
1029 y0 = sVloc[2];
1030 y0a = -1.;
1031 i1 = (sNpy << 1) + 1;
1032 for (i = 1; i <= i1; ++i) {
1033 if (y0 > sYanode[i - 1] && y0 <= sYanode[i]) {
1034 ass_i1 = TMath::Abs(sYanode[i] - y0);
1035 ass_i = TMath::Abs(sYanode[i - 1] - y0);
1036 if (ass_i1 <= ass_i) y0a = sYanode[i];
1037 if (ass_i1 > ass_i) y0a = sYanode[i - 1];
1038 return;
1039 }
1040 }
1041}
1042
1043//_____________________________________________________________________________
1044void AliRICH::GetChargeMip(Float_t &qtot)
1045{
1046 //
1047 // Calculates the charge deposited by a MIP
1048 //
1049
1050 AliMC* pMC = AliMC::GetMC();
1051
1052 Int_t i1;
1053
1054 // Local variables
1055 Float_t ranf[2];
1056 Int_t i;
1057 Int_t nel;
1058
1059 // NUMBER OF ELECTRONS
1060
1061 nel = (Int_t) (sDetot * 1e9 / 26.);
1062 qtot = 0.;
1063 if (!nel) nel = 1;
1064 i1 = nel;
1065 for (i = 1; i <= i1; ++i) {
1066 pMC->Rndm(ranf, 1);
1067 qtot -= fChslope * TMath::Log(ranf[0]);
1068 }
1069}
1070
1071//_____________________________________________________________________________
1072void AliRICH::GetCharge(Float_t &qtot)
1073{
1074 //
1075 // Charge deposited
1076 //
1077
1078 AliMC* pMC = AliMC::GetMC();
1079
1080 Float_t ranf[1];
1081
1082 pMC->Rndm(ranf, 1);
1083 qtot = -fChslope * TMath::Log(ranf[0]);
1084}
1085
1086//_____________________________________________________________________________
1087void AliRICH::FeedBack(Float_t *source, Float_t qtot)
1088{
1089 //
1090 // Generate FeedBack photons
1091 //
1092
1093 AliMC* pMC = AliMC::GetMC();
1094 TGeant3 *geant3 = (TGeant3*) pMC;
1095
1096 Int_t i1, j;
1097 Float_t r1, r2;
1098
1099 // Local variables
1100 Float_t cthf, ranf[2], phif, enfp = 0, sthf;
1101 Int_t i, ifeed;
1102 Float_t e1[3], e2[3], e3[3];
1103 Float_t vmod, uswop;
1104 Float_t fp, random;
1105 Float_t dir[3], phi;
1106 Int_t nfp;
1107 Float_t pol[3], supwght;
1108
1109 // DETERMINE NUMBER OF FEEDBACK PHOTONS
1110
1111 // Function Body
1112 fp = fAlphaFeed * qtot;
1113 nfp = gRandom->Poisson(fp);
1114
1115 // GENERATE PHOTONS
1116
1117 geant3->Gckin2()->ngphot = 0;
1118 i1 = nfp;
1119 for (i = 1; i <= i1; ++i) {
1120
1121 // DIRECTION
1122 pMC->Rndm(ranf, 2);
1123 cthf = ranf[0] * 2. - 1.;
1124 if (cthf < 0.) continue;
1125 sthf = TMath::Sqrt((1. - cthf) * (cthf + 1.));
1126 phif = ranf[1] * 2 * TMath::Pi();
1127 ++geant3->Gckin2()->ngphot;
1128 if (geant3->Gckin2()->ngphot > 800) {
1129 printf("ATTENTION NGPHOT ! %d %f %d\n",
1130 geant3->Gckin2()->ngphot,fp,nfp);
1131 break;
1132 }
1133
1134 // POSITION
1135 for(j=0;j<3;j++)
1136 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][j]=source[j];
1137
1138 // ENERGY
1139 pMC->Rndm(&random, 1);
1140 if (random <= .57) {
1141 enfp = 7.5e-9;
1142 } else if (random > .57 && random <= .7) {
1143 enfp = 6.4e-9;
1144 } else if (random > .7) {
1145 enfp = 7.9e-9;
1146 }
1147 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][6] = enfp;
1148 dir[0] = sthf * TMath::Sin(phif);
1149 dir[1] = cthf;
1150 dir[2] = sthf * TMath::Cos(phif);
1151 pMC->Gdtom(dir, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][3], 2);
1152
1153 // POLARISATION
1154 e1[0] = 0.;
1155 e1[1] = -dir[2];
1156 e1[2] = dir[1];
1157
1158 e2[0] = -dir[1];
1159 e2[1] = dir[0];
1160 e2[2] = 0.;
1161
1162 e3[0] = dir[1];
1163 e3[1] = 0.;
1164 e3[2] = -dir[0];
1165
1166 vmod=0;
1167 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
1168 if (!vmod) for(j=0;j<3;j++) {
1169 uswop=e1[j];
1170 e1[j]=e3[j];
1171 e3[j]=uswop;
1172 }
1173 vmod=0;
1174 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
1175 if (!vmod) for(j=0;j<3;j++) {
1176 uswop=e2[j];
1177 e2[j]=e3[j];
1178 e3[j]=uswop;
1179 }
1180
1181 vmod=0;
1182 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
1183 vmod=TMath::Sqrt(1/vmod);
1184 for(j=0;j<3;j++) e1[j]*=vmod;
1185
1186 vmod=0;
1187 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
1188 vmod=TMath::Sqrt(1/vmod);
1189 for(j=0;j<3;j++) e2[j]*=vmod;
1190
1191 pMC->Rndm(ranf, 1);
1192 phi = ranf[0] * 2 * TMath::Pi();
1193 r1 = TMath::Sin(phi);
1194 r2 = TMath::Cos(phi);
1195 for(j=0;j<3;j++) pol[j]=e1[j]*r1+e2[j]*r2;
1196 pMC->Gdtom(pol, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][7], 2);
1197
1198 // TIME OF FLIGHT
1199 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][10] = 0.;
1200
1201 // PUT PHOTON ON THE STACK AND LABEL IT AS FEEDBACK (51,52)
1202 supwght = geant3->Gctrak()->upwght;
1203 r1 = geant3->Gctrak()->upwght / 100.;
1204 ifeed = Int_t(r1+0.5);
1205 ++sNfeed;
1206 if (geant3->Gckine()->ipart == 50 && ifeed != 52) {
1207 geant3->Gctrak()->upwght = 5100.;
1208 } else {
1209 geant3->Gctrak()->upwght = 5200.;
1210 }
1211 geant3->Gskpho(geant3->Gckin2()->ngphot);
1212 geant3->Gctrak()->upwght = supwght;
1213
1214 }
1215 geant3->Gckin2()->ngphot = 0;
1216}
1217
1218//_____________________________________________________________________________
1219Float_t AliRICH::FMathieson(Float_t lambda1, Float_t lambda2)
1220{
1221 //
1222 // Mathieson charge distribution
1223 //
1224 // CALCULATES INTEGRAL OF GATTI/MATHIESON CHARGE DISTRIBUTION
1225 // FOR CPC AND CSC CHAMBERS
1226 // INTEGRATION LIMITS LAMBDA1,LAMBDA2= X/D WHERE D IS THE ANODE CATHODE
1227 // SEPARATION
1228 // K3: GEOMETRY PARAMETER
1229 //
1230 Float_t u1, u2;
1231 // const Float_t k1=0.2828;
1232 const Float_t k2=0.962;
1233 const Float_t k3=0.6;
1234 const Float_t k4=0.379;
1235 const Float_t sqrtk3=TMath::Sqrt(k3);
1236
1237
1238 u1 = tanh(lambda1 * k2) * sqrtk3;
1239 u2 = tanh(lambda2 * k2) * sqrtk3;
1240
1241 return (atan(u2) - atan(u1)) * 2 * k4;
1242
1243}
1244
1245
1246//_____________________________________________________________________________
1247void AliRICH::UpdateMipHit(Float_t *hits)
1248{
1249 //
1250 // Update some parameters when the current mip hits the detector.
1251 //
1252 TClonesArray &lhits = *fMips;
1253 AliRICHmip *lmip = (AliRICHmip *) lhits.AddrAt(fNmips-1);
1254 lmip->SetX ( hits[1]);
1255 lmip->SetY ( hits[2]);
1256 lmip->SetModule((int) hits[3]);
1257 lmip->SetTheta ( hits[4]);
1258 lmip->SetPhi ( hits[5]);
1259}
1260
1261//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
1262
1263//_____________________________________________________________________________
1264void AliRICH::MakeBranch(Option_t *option)
1265{
1266 //
1267 // Create a new branch in the current Root Tree.
1268 // The branch of fHits is automatically split
1269 //
1270 AliDetector::MakeBranch(option);
1271
1272 Int_t buffersize = 4000;
1273 if (gAlice->TreeH()) {
1274 if(fMips){
1275 gAlice->TreeH()->Branch("RICHmip",&fMips, buffersize);
1276 printf("Making Branch %s for mips\n","RICHmip");
1277 }
1278 if(fCkovs){
1279 gAlice->TreeH()->Branch("RICHckov",&fCkovs, buffersize);
1280 printf("Making Branch %s for ckovs\n","RICHckov");
1281 }
1282 if(fPadhits){
1283 gAlice->TreeH()->Branch("RICHpadhit",&fPadhits, buffersize);
1284 printf("Making Branch %s for padhits\n","RICHpadhit");
1285 }
1286 }
1287}
1288
1289//_____________________________________________________________________________
1290void AliRICH::SetTreeAddress()
1291{
1292 //
1293 // Set branch address for the Hits and Digits Tree.
1294 //
1295 AliDetector::SetTreeAddress();
1296 TBranch *branch;
1297 TTree *treeH = gAlice->TreeH();
1298 if (treeH) {
1299 if (fMips) {
1300 branch = treeH->GetBranch("RICHmip");
1301 if (branch) branch->SetAddress(&fMips);
1302 }
1303 if (fCkovs) {
1304 branch = treeH->GetBranch("RICHckov");
1305 if (branch) branch->SetAddress(&fCkovs);
1306 }
1307 if (fPadhits) {
1308 branch = treeH->GetBranch("RICHpadhit");
1309 if (branch) branch->SetAddress(&fPadhits);
1310 }
1311 }
1312}
1313
1314//_____________________________________________________________________________
1315void AliRICH::ResetHits()
1316{
1317 //
1318 // Reset number of mips and the mips array for RICH
1319 //
1320 AliDetector::ResetHits();
1321 fNmips = 0;
1322 if (fMips) fMips->Clear();
1323 // Reset number of Ckovs and the Ckovs array for RICH
1324 fNckovs = 0;
1325 if (fCkovs) fCkovs->Clear();
1326 // Reset number of Padhits and the Padhits array for RICH
1327 fNpadhits = 0;
1328 if (fPadhits) fPadhits->Clear();
1329}
1330
1331//_____________________________________________________________________________
1332Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1333{
1334 //
1335 // Distance from mouse RICH on the display
1336 // Dummy routine
1337 //
1338 return 9999;
1339}
1340
1341//_____________________________________________________________________________
1342void AliRICH::PreTrack()
1343{
1344 //
1345 // Called before transporting a track
1346 //
1347 sDetot=0;
1348 sNsecon=0;
1349 sNpads=0;
1350 sNphoton=0;
1351 sNfeed=0;
1352 sNfeedd=0;
1353 sNdir=0;
1354}
1355
1356//_____________________________________________________________________________
1357void AliRICH::Streamer(TBuffer &R__b)
1358{
1359 //
1360 // Stream an object of class AliRICH.
1361 //
1362 if (R__b.IsReading()) {
1363 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1364 AliDetector::Streamer(R__b);
1365 R__b >> fNmips;
1366 R__b >> fNckovs;
1367 R__b >> fNpadhits;
1368 R__b >> fMips; //diff
1369 R__b >> fCkovs; //diff
1370 R__b >> fPadhits; //diff
1371 R__b >> fChslope;
1372 R__b >> fAlphaFeed;
1373 R__b >> fSxcharge;
1374 R__b >> fIritri;
1375 } else {
1376 R__b.WriteVersion(AliRICH::IsA());
1377 AliDetector::Streamer(R__b);
1378 R__b << fNmips;
1379 R__b << fNckovs;
1380 R__b << fNpadhits;
1381 R__b << fMips; //diff
1382 R__b << fCkovs; //diff
1383 R__b << fPadhits; //diff
1384 R__b << fChslope;
1385 R__b << fAlphaFeed;
1386 R__b << fSxcharge;
1387 R__b << fIritri;
1388 }
1389}
1390
1391ClassImp(AliRICHv1)
1392
1393///////////////////////////////////////////////////////////////////////////////
1394// //
1395// Ring Imaging Cherenkov //
1396// This class contains version 1 of the Ring Imaging Cherenkov //
1397// //
1398//Begin_Html
1399/*
1400<img src="gif/AliRICHv1Class.gif">
1401*/
1402//End_Html
1403// //
1404///////////////////////////////////////////////////////////////////////////////
1405
1406//_____________________________________________________________________________
1407AliRICHv1::AliRICHv1() : AliRICH()
1408{
1409 //
1410 // Default Constructor RICH for version 1
1411 //
1412}
1413
1414//_____________________________________________________________________________
1415AliRICHv1::AliRICHv1(const char *name, const char *title)
1416 : AliRICH(name,title)
1417{
1418 //
1419 // Standard constructor for RICH version 1
1420 //
1421}
1422
1423//_____________________________________________________________________________
1424AliRICHv1::~AliRICHv1()
1425{
1426 //
1427 // Standard destructor for RICH version 1
1428 //
1429}
1430
1431//_____________________________________________________________________________
1432void AliRICHv1::CreateGeometry()
1433{
1434 //
1435 // Create the geometry for RICH version 1
1436 //
1437 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1438 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
1439 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1440 //
1441 //Begin_Html
1442 /*
1443 <img src="gif/AliRICHv1.gif">
1444 */
1445 //End_Html
1446 //Begin_Html
1447 /*
1448 <img src="gif/AliRICHv1Tree.gif">
1449 */
1450 //End_Html
1451
1452 AliMC* pMC = AliMC::GetMC();
1453
1454 Int_t *idtmed = gAlice->Idtmed();
1455
1456 Int_t i;
1457 Float_t zs;
1458 Int_t idrotm[1099];
1459 Float_t par[3];
1460
1461 // --- Define the RICH detector
1462 // External aluminium box
1463 par[0] = 71.1;
1464 par[1] = 11.5;
1465 par[2] = 73.15;
1466 pMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
1467
1468 // Sensitive part of the whole RICH
1469 par[0] = 64.8;
1470 par[1] = 11.5;
1471 par[2] = 66.55;
1472 pMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
1473
1474 // Honeycomb
1475 par[0] = 63.1;
1476 par[1] = .188;
1477 par[2] = 66.55;
1478 pMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
1479
1480 // Aluminium sheet
1481 par[0] = 63.1;
1482 par[1] = .025;
1483 par[2] = 66.55;
1484 pMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
1485
1486 // Quartz
1487 par[0] = 63.1;
1488 par[1] = .25;
1489 par[2] = 65.5;
1490 pMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
1491
1492 // Spacers (cylinders)
1493 par[0] = 0.;
1494 par[1] = .5;
1495 par[2] = .5;
1496 pMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
1497
1498 // Opaque quartz
1499 par[0] = 61.95;
1500 par[1] = .2;
1501 par[2] = 66.5;
1502 pMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
1503
1504 // Frame of opaque quartz
1505 par[0] = 20.65;
1506 par[1] = .5;
1507 par[2] = 66.5;
1508 pMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
1509
1510 // Little bar of opaque quartz
1511 par[0] = 63.1;
1512 par[1] = .25;
1513 par[2] = .275;
1514 pMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
1515
1516 // Freon
1517 par[0] = 20.15;
1518 par[1] = .5;
1519 par[2] = 65.5;
1520 pMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
1521
1522 // Methane
1523 par[0] = 64.8;
1524 par[1] = 5.;
1525 par[2] = 64.8;
1526 pMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
1527
1528 // Methane gap
1529 par[0] = 64.8;
1530 par[1] = .2;
1531 par[2] = 64.8;
1532 pMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1533
1534 // CsI photocathode
1535 par[0] = 64.8;
1536 par[1] = .25;
1537 par[2] = 64.8;
1538 pMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1539
1540 // Anode grid
1541 par[0] = 0.;
1542 par[1] = .0025;
1543 par[2] = 20.;
1544 pMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1545
1546 // --- Places the detectors defined with GSVOLU
1547 // Place material inside RICH
1548 pMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
1549
1550 pMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
1551 pMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
1552 pMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
1553 pMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
1554
1555 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1556
1557 for (i = 1; i <= 9; ++i) {
1558 zs = (5 - i) * 14.4;
1559 pMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
1560 }
1561 for (i = 10; i <= 18; ++i) {
1562 zs = (14 - i) * 14.4;
1563 pMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
1564 }
1565
1566 pMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
1567 pMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
1568 pMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
1569 pMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
1570 pMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
1571 pMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
1572 pMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
1573 pMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
1574 pMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1575 pMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
1576
1577 // Place RICH inside ALICE apparatus
1578
1579 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1580 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1581 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1582 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1583 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1584 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1585 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1586
1587 pMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
1588 pMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
1589 pMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
1590 pMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
1591 pMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
1592 pMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
1593 pMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
1594
1595}
1596
1597//_____________________________________________________________________________
860ec329 1598void AliRICHv1::DrawModule()
fe4da5cc 1599{
1600 //
1601 // Draw a shaded view of the Ring Imaging Cherenkov
1602 //
1603
1604 AliMC* pMC = AliMC::GetMC();
1605 TGeant3 *geant3 = (TGeant3*) pMC;
1606
1607 // Set everything unseen
1608 pMC->Gsatt("*", "seen", -1);
1609 //
1610 // Set ALIC mother transparent
1611 pMC->Gsatt("ALIC","SEEN",0);
1612 //
1613 // Set the volumes visible
1614
1615 pMC->Gsatt("RICH","seen",0);
1616 pMC->Gsatt("SRIC","seen",0);
1617 pMC->Gsatt("HONE","seen",1);
1618 pMC->Gsatt("ALUM","seen",1);
1619 pMC->Gsatt("QUAR","seen",1);
1620 pMC->Gsatt("SPAC","seen",1);
1621 pMC->Gsatt("OQUA","seen",1);
1622 pMC->Gsatt("OQUF","seen",1);
1623 pMC->Gsatt("BARR","seen",1);
1624 pMC->Gsatt("FREO","seen",1);
1625 pMC->Gsatt("META","seen",1);
1626 pMC->Gsatt("GAP ","seen",1);
1627 pMC->Gsatt("CSI ","seen",1);
1628 pMC->Gsatt("GRID","seen",1);
1629
1630 //
1631 geant3->Gdopt("hide", "on");
1632 geant3->Gdopt("shad", "on");
1633 geant3->Gsatt("*", "fill", 7);
1634 geant3->SetClipBox(".");
1635 geant3->Gdopt("hide", "on");
1636 geant3->Gdopt("shad", "on");
1637 geant3->Gsatt("*", "fill", 7);
1638 geant3->SetClipBox(".");
1639 // geant3->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
1640 geant3->DefaultRange();
1641 pMC->Gdraw("alic", 60, 50, 0, 10, 0, .03, .03);
1642 geant3->Gdhead(1111, "Ring Imaging Cherenkov version 1");
1643 geant3->Gdman(16, 6, "MAN");
1644 geant3->Gdopt("hide", "off");
1645}
1646
1647//_____________________________________________________________________________
1648void AliRICHv1::CreateMaterials()
1649{
1650 //
1651 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1652 // ORIGIN : NICK VAN EIJNDHOVEN
1653 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1654 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
1655 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1656 //
1657 Int_t ISXFLD = gAlice->Field()->Integ();
1658 Float_t SXMGMX = gAlice->Field()->Max();
1659
1660 Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
1661 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1662 Float_t rindex_quarz[14] = { 1.528309,1.533333,
1663 1.538243,1.544223,1.550568,1.55777,
1664 1.565463,1.574765,1.584831,1.597027,
1665 1.611858,1.6277,1.6472,1.6724 };
1666 Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1667 Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1668 Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1669 Float_t absco_freon[14] = { 179.0987,179.0987,
1670 179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
1671 Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
1672 6.324,4.483,1.6,.323,.073,0. };
1673 Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1674 1e-5,1e-5,1e-5,1e-5,1e-5 };
1675 Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1676 1e-4,1e-4,1e-4,1e-4 };
1677 Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1678 1e6,1e6,1e6 };
1679 Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1680 1e-4,1e-4,1e-4,1e-4 };
1681 Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1682 Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
1683 .2101,.2554,.293,.376,.3861,.418 };
1684 Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1685
1686 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1687 zmet[2], zqua[2];
1688 Int_t nlmatfre;
1689 Float_t densquao;
1690 Int_t nlmatmet, nlmatqua;
1691 Float_t wmatquao[2], rindex_freon[14];
1692 Int_t i;
1693 Float_t aquao[2], epsil, stmin, zquao[2];
1694 Int_t nlmatquao;
1695 Float_t radlal, densal, tmaxfd, deemax, stemax;
1696 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1697
1698 Int_t *idtmed = gAlice->Idtmed();
1699
1700 AliMC* pMC = AliMC::GetMC();
1701 TGeant3 *geant3 = (TGeant3*) pMC;
1702
1703 // --- Photon energy (GeV)
1704 // --- Refraction indexes
1705 for (i = 0; i < 14; ++i) {
1706 rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
1707 }
1708 // need to be changed
1709
1710 // --- Absorbtion lenghts (in cm)
1711 // DATA ABSCO_QUARZ /
1712 // & 5 * 1000000.,
1713 // & 29.85, 7.34, 4.134, 1.273, 0.722,
1714 // & 0.365, 0.365, 0.365, 0. /
1715 // need to be changed
1716
1717 // --- Detection efficiencies (quantum efficiency for CsI)
1718 // --- Define parameters for honeycomb.
1719 // Used carbon of equivalent rad. lenght
1720
1721 ahon = 12.01;
1722 zhon = 6.;
1723 denshon = 2.265;
1724 radlhon = 18.8;
1725
1726 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1727
1728 aqua[0] = 28.09;
1729 aqua[1] = 16.;
1730 zqua[0] = 14.;
1731 zqua[1] = 8.;
1732 densqua = 2.64;
1733 nlmatqua = -2;
1734 wmatqua[0] = 1.;
1735 wmatqua[1] = 2.;
1736
1737 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1738
1739 aquao[0] = 28.09;
1740 aquao[1] = 16.;
1741 zquao[0] = 14.;
1742 zquao[1] = 8.;
1743 densquao = 2.64;
1744 nlmatquao = -2;
1745 wmatquao[0] = 1.;
1746 wmatquao[1] = 2.;
1747
1748 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1749
1750 afre[0] = 12.;
1751 afre[1] = 19.;
1752 zfre[0] = 6.;
1753 zfre[1] = 9.;
1754 densfre = 1.7;
1755 nlmatfre = -2;
1756 wmatfre[0] = 6.;
1757 wmatfre[1] = 14.;
1758
1759 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1760
1761 amet[0] = 12.01;
1762 amet[1] = 1.;
1763 zmet[0] = 6.;
1764 zmet[1] = 1.;
1765 densmet = 7.17e-4;
1766 nlmatmet = -2;
1767 wmatmet[0] = 1.;
1768 wmatmet[1] = 4.;
1769
1770 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1771
1772 agri = 63.54;
1773 zgri = 29.;
1774 densgri = 8.96;
1775 radlgri = 1.43;
1776
1777 // --- Parameters to include in GSMATE related to aluminium sheet
1778
1779 aal = 26.98;
1780 zal = 13.;
1781 densal = 2.7;
1782 radlal = 8.9;
1783
1784 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1785 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1786 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1787 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1788 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1789 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1790 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1791 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1792 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1793 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1794
1795 tmaxfd = -10.;
1796 stemax = -.1;
1797 deemax = -.2;
1798 epsil = .001;
1799 stmin = -.001;
1800
1801 AliMedium(1001, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1802 AliMedium(1002, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1803 AliMedium(1003, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1804 AliMedium(1004, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1805 AliMedium(1005, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1806 AliMedium(1006, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
1807 AliMedium(1007, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1808 AliMedium(1008, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1809 AliMedium(1009, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
1810 AliMedium(1010, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1811
1812
1813 // Switch on delta-ray production in the methane and freon gaps
1814
1815 pMC->Gstpar(idtmed[1002], "LOSS", 1.);
1816 pMC->Gstpar(idtmed[1003], "LOSS", 1.);
1817 pMC->Gstpar(idtmed[1004], "LOSS", 1.);
1818 pMC->Gstpar(idtmed[1008], "LOSS", 1.);
1819 pMC->Gstpar(idtmed[1005], "LOSS", 1.);
1820 pMC->Gstpar(idtmed[1002], "HADR", 1.);
1821 pMC->Gstpar(idtmed[1003], "HADR", 1.);
1822 pMC->Gstpar(idtmed[1004], "HADR", 1.);
1823 pMC->Gstpar(idtmed[1008], "HADR", 1.);
1824 pMC->Gstpar(idtmed[1005], "HADR", 1.);
1825 pMC->Gstpar(idtmed[1002], "DCAY", 1.);
1826 pMC->Gstpar(idtmed[1003], "DCAY", 1.);
1827 pMC->Gstpar(idtmed[1004], "DCAY", 1.);
1828 pMC->Gstpar(idtmed[1008], "DCAY", 1.);
1829 pMC->Gstpar(idtmed[1005], "DCAY", 1.);
1830 geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
1831 geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
1832 geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
1833 geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
1834 geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
1835 geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
1836 geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
1837 geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
1838 geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
1839 geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
1840}
1841
1842ClassImp(AliRICHhit)
1843
1844//_____________________________________________________________________________
1845AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol,
1846 Float_t *hits, Int_t fNpadhits) :
1847 AliHit(shunt,track)
1848{
1849 //
1850 // Standard constructor for RICH hit
1851 //
1852 Int_t i;
1853 for (i=0;i<2;i++) fVolume[i] = vol[i];
1854 fTrack = (int) track;
1855 //Hit information
1856 fPart = (int) hits[ 1];
1857 //AliHit information, position of the hit
1858 fX = hits[ 2];
1859 fY = hits[ 3];
1860 //Pad information
1861 fFirstpad = fNpadhits;
1862 fLastpad = fNpadhits-1;
1863
1864 //Hit information
1865 fModule = (int) hits[ 4];
1866 fTheta = hits[ 5];
1867 fArrivaltime= hits[ 6];
1868 fFeed = (int) hits[ 7];
1869}
1870
1871ClassImp(AliRICHmip)
1872
1873//_____________________________________________________________________________
1874AliRICHmip::AliRICHmip(Int_t shunt, Int_t track, Int_t *vol,
1875 Float_t *hits, Int_t fNckovs, Int_t fNpadhits) :
1876 AliRICHhit(shunt,track,vol,hits,fNpadhits)
1877{
1878 //
1879 // Standard constructor for RICH Mip hit
1880 //
1881 fPhi = hits[ 8];
1882 fPs = hits[ 9];
1883 fQ = hits[10];
1884 fZ = hits[11];
1885 //Ckov information
1886 if ((int) hits[12]){
1887 fFirstCkov = fNckovs;
1888 fLastCkov = fNckovs-1;
1889 } else {
1890 fFirstCkov = -1;
1891 fLastCkov = -1;
1892 }
1893}
1894
1895ClassImp(AliRICHckov)
1896
1897//_____________________________________________________________________________
1898AliRICHckov::AliRICHckov(Int_t shunt, Int_t track, Int_t *vol,
1899 Float_t *hits, Int_t fNmips, Int_t fNpadhits) :
1900 AliRICHhit(shunt,track,vol,hits,fNpadhits)
1901{
1902 //
1903 // Standard creator for RICH Cherenkov information
1904 //
1905 fEnergy = hits[8];
1906 fStop = (int) hits[9];
1907
1908 //Parent info
1909 fParent = fNmips;
1910}
1911
1912ClassImp(AliRICHpadhit)
1913
1914
1915//_____________________________________________________________________________
1916AliRICHpadhit::AliRICHpadhit(Int_t shunt, Int_t track, Int_t *vol,
1917 Float_t *hits, Int_t fNmips, Int_t fNckovs):
1918 AliHit(shunt,track)
1919{
1920 //
1921 // Standard constructor for RICH pad hit
1922 //
1923 Int_t i;
1924 for (i=0;i<2;i++) fVolume[i] = vol[i];
1925
1926 // Hit information
1927 fX = (int) hits[ 2];
1928 fY = (int) hits[ 3];
1929 fModule = (int) hits[ 4];
1930 fParentMip = fNmips;
1931 fParentCkov = fNckovs;
1932 fProd = (int) hits[ 5];
1933 fCharge = hits[ 6];
1934 fZ = -999.0; // Not implemented
1935}
1936
1937