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