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