1 ///////////////////////////////////////////////////////////////////////////////
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 //
10 <img src="picts/AliRICHClass.gif">
15 ///////////////////////////////////////////////////////////////////////////////
29 //_____________________________________________________________________________
33 // Default constructor for RICH
50 //_____________________________________________________________________________
51 AliRICH::AliRICH(const char *name, const char *title)
52 : AliDetector(name,title)
55 // Standard constructor for RICH
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);
68 // Set parameters to default value
79 //_____________________________________________________________________________
83 // Destructor for RICH
87 fMips->Delete(); delete fMips;
88 fCkovs->Delete(); delete fCkovs;
89 fPadhits->Delete(); delete fPadhits;
92 //_____________________________________________________________________________
93 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
98 switch ((int) hits[0]) {
103 TClonesArray &lhits = *fHits;
104 new(lhits[fNhits++]) AliHit(fIshunt,track);
105 AliHit *lhit = (AliHit*) fHits->AddrAt(fNhits-1);
113 AddMipHit(track,vol,hits);
118 AddCkovHit(track,vol,hits);
123 AddPadHit(track,vol,hits);
130 printf("Error: AliRICH::AddHit flag %d not defined./n",(int) hits[0]);
135 //_____________________________________________________________________________
136 void AliRICH::AddMipHit(Int_t track, Int_t *vol, Float_t *hits)
138 // Adds a mip hit in the RICH.
139 TClonesArray &lhits = *fMips;
140 new(lhits[fNmips++]) AliRICHmip(fIshunt,track,vol,hits,
144 //_____________________________________________________________________________
145 void AliRICH::AddCkovHit(Int_t track, Int_t *vol, Float_t *hits)
148 // Adds a cerenkov hit in the RICH.
150 TClonesArray &lhits = *fCkovs;
151 AliRICHmip *lmip = (AliRICHmip*) fMips->AddrAt(fNmips-1);
153 // If this ckov come from a mip update the mip lastckov entry.
155 if (lmip->GetZ() != -999.0) {
156 fmipslocal = fNmips-1;
157 lmip->SetLastCkov(fNckovs);
159 new(lhits[fNckovs++]) AliRICHckov(fIshunt,track,vol,hits,
160 fmipslocal,fNpadhits);
163 //_____________________________________________________________________________
164 void AliRICH::AddPadHit(Int_t track, Int_t *vol, Float_t *hits)
167 // Adds pad hits in the RICH
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);
175 ((AliRICHmip *) fMips->AddrAt(fNmips-1))->SetLastpad(fNpadhits);
176 new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,fNmips-1,-1);
180 //_____________________________________________________________________________
181 void AliRICH::BuildGeometry()
184 // Builds a TNode geometry for event display
188 const int kColorRICH = kGreen;
190 Top=gAlice->GetGeometry()->GetNode("alice");
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);
201 Node = new TNode("RICH1","RICH1","S_RICH",0,471.8999,165.2599,"rot993");
202 Node->SetLineColor(kColorRICH);
205 Node = new TNode("RICH2","RICH2","S_RICH",171,470,0,"rot994");
206 Node->SetLineColor(kColorRICH);
209 Node = new TNode("RICH3","RICH3","S_RICH",0,500,0,"rot995");
210 Node->SetLineColor(kColorRICH);
213 Node = new TNode("RICH4","RICH4","S_RICH",-171,470,0,"rot996");
214 Node->SetLineColor(kColorRICH);
217 Node = new TNode("RICH5","RICH5","S_RICH",161.3999,443.3999,-165.3,"rot997");
218 Node->SetLineColor(kColorRICH);
221 Node = new TNode("RICH6","RICH6","S_RICH",0,471.8999,-165.3,"rot998");
222 Node->SetLineColor(kColorRICH);
225 Node = new TNode("RICH7","RICH7","S_RICH",-161.399,443.3999,-165.3,"rot999");
226 Node->SetLineColor(kColorRICH);
230 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
232 // Section for Rich Step
235 static Float_t sVloc[3];
236 static Float_t sVectIn[3];
237 static Float_t sDetot;
238 static Float_t sYanode[MAXPH];
239 static Float_t sXpad[MAXPH];
240 static Float_t sYpad[MAXPH];
247 static Float_t sDpad;
250 static Int_t sNsecon;
251 //static Float_t sQint[2];
256 static Int_t sIpx[MAXSEC];
257 static Int_t sIpy[MAXSEC];
258 static Int_t sIqpad[MAXSEC];
259 static Int_t sNphlink[MAXSEC];
260 static Int_t sNphoton;
262 static Int_t sNfeed, sNfeedd, sNdir;
264 static Float_t sPparent;
265 static Float_t sThincParent;
266 static Int_t sIloss[MAXPH];
267 static Int_t sIprod[MAXPH];
268 static Float_t sXphit[MAXPH];
269 static Float_t sYphit[MAXPH];
271 static Float_t sSycharge;
272 static Float_t sXsox, sYsox, sZsox;
274 static Int_t sMckov[MAXPH];
275 static Int_t idpartgx;
276 static Float_t phisx;
281 static Float_t thetasx;
284 //_____________________________________________________________________________
288 // Initialise RICH after that it has been built
290 const Float_t sconv=2*TMath::Sqrt(2*TMath::Log(2));
291 const Float_t yK3=1.20;
304 for(i=0;i<sNpx;i++) {
308 for(i=0;i<2*sNpy+1;i++) sYanode[i]=ansp/2+i*ansp;
311 sSycharge=4*TMath::ATanH(1/TMath::Sqrt(2*yK3))/TMath::Pi()
312 /(1-0.5*TMath::Sqrt(yK3))/sDpad/sconv;
317 //___________________________________________________________________________
318 void AliRICH::StepManager()
321 // Called at every step in the RICH
324 TGeant3 *geant3 = (TGeant3*) gMC;
326 const Float_t xshift[3] = { 41.3, 0, -41.3 };
327 const Int_t nrooth = 25;
329 static Int_t ixold=-1, iyold=-1;
331 // System generated locals
336 Float_t ranf[2], rrhh[nrooth], phiangle, cost, vmod;
351 //Float_t xtrig[200], ytrig[200];
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
364 // IPROD = 1 PRODUCED IN FREON
365 // IPROD = 2 PRODUCED IN QUARZ
367 // new (changed NROOTH from 10 to 25!!!!!!!!!!!!!)
370 Int_t *idtmed = fIdtmed->GetArray()-999;
372 //--------------------------------------------------------------------------
376 if (geant3->Gckine()->charge) {
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) {
386 gMC->Gmtod(geant3->Gckin3()->gpos[i-1], vxloc, 1);
387 gMC->Gmtod(geant3->Gcking()->gkin[i-1], dir, 2);
388 if (geant3->Gcking()->gkin[i-1][4] == 8. ||
389 geant3->Gcking()->gkin[i-1][4] == 9.) {
391 // Computing 2nd power
393 // Computing 2nd power
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];
404 if ((geant3->Gctmed()->numed == idtmed[1006-1] &&
405 geant3->Gctrak()->inwvol == 2) ||
406 geant3->Gctrak()->istop) {
408 printf("NOT TRIGGERED\n");
416 geant3->Gctrak()->istory = 0;
417 geant3->Gctrak()->upwght = 0.;
418 geant3->Gcflag()->ieotri = 1;
423 printf("TRIGGERED %d\n",nmult);
427 // MIP inside Methane
428 if (geant3->Gctmed()->numed == idtmed[1009-1]) {
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 ?
435 gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
436 gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
437 phiangle = TMath::ATan2(dir[2], dir[0]);
438 if (phiangle < 0.) phiangle += 2*TMath::Pi();
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];
444 ihitrak = gAlice->CurrentTrack();
446 // flag to say this is update
447 rrhh[1] = sVloc[0] + sDlx;
449 rrhh[2] = sVloc[2] + sDly;
451 rrhh[3] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
453 // Computing 2nd power
455 // Computing 2nd power
457 rrhh[4] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
460 // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
461 AddHit(ihitrak,irivol,rrhh);
464 // Record particle properties
465 // If particle produced already Cerenkov Photons (istory=1)
467 // update the impact point only
468 if (geant3->Gctrak()->istory != 2) {
469 if (!geant3->Gctrak()->istory) {
472 // Is this a primary particle ?
474 //if (geant3->Gctrak()->upwght) iprimx = 0;
476 // Where did it come from ?
477 sXsox = geant3->Gckine()->vert[0];
478 sYsox = geant3->Gckine()->vert[1];
479 sZsox = geant3->Gckine()->vert[2];
482 psx = geant3->Gctrak()->vect[6];
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;
493 // Direction of incidence and where did it hit ?
494 gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
495 gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
496 // Computing 2nd power
498 // Computing 2nd power
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];
508 for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
509 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
511 ihitrak = gAlice->CurrentTrack();
513 // Flag to say this is MIP
514 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
517 rrhh[4] = (Float_t) nmodsx;
520 rrhh[6] = geant3->Gctrak()->tofg;
522 rrhh[7] = (Float_t) idpartgx;
526 // charge of current particle in electron charge unit;
527 rrhh[10] = geant3->Gckine()->charge;
529 // Zo of ckov generation
532 AddHit(ihitrak, irivol, rrhh);
535 // Earmark track as being recorded in methane gap
537 geant3->Gctrak()->istory = 2;
541 // Signal generation in methane gap
542 gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
543 gMC->Gmtod(geant3->Gckine()->vert, vxloc, 1);
544 ix = (Int_t) ((sVloc[0] + sDlx) / sDxp);
545 iy = (Int_t) ((sVloc[2] + sDly) / sDyp);
547 // Is this the first step?
548 if (ixold == -1 && iyold == -1) {
551 for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
555 if (geant3->Gctrak()->inwvol == 2 || geant3->Gctrak()->istop) {
556 sDetot += geant3->Gctrak()->destep;
557 if (sDetot > 0.) RichIntegration();
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;
570 sDetot += geant3->Gctrak()->destep;
575 // End charged particles treatment
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) {
583 // If it is a primary, save all generated photons
584 i1 = geant3->Gckin2()->ngphot;
585 for (i = 1; i <= i1; ++i) {
587 if (sNphoton > MAXPH) {
589 printf("ATTENTION NPHOTON %d\n",sNphoton);
595 if (geant3->Gctmed()->numed == idtmed[1003-1]) medprod = 2;
597 // Production angle and energy
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];
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;
612 // geant3->Gskpho(i);
613 sMckov[sNphoton - 1] = gAlice->CurrentTrack()+i;
614 geant3->Gctrak()->upwght = stwght;
617 stwght = geant3->Gctrak()->upwght;
618 geant3->Gctrak()->upwght = 0.;
619 // geant3->Gskpho(0);
620 geant3->Gctrak()->upwght = stwght;
623 // Particle did not yet pass the methane gap
624 if (geant3->Gctrak()->istory == 0) {
625 geant3->Gctrak()->istory = 1;
627 // Is this a primary particle ?
629 //if (geant3->Gctrak()->upwght) iprimx = 0;
631 // Where did it come from ?
632 sXsox = geant3->Gckine()->vert[0];
633 sYsox = geant3->Gckine()->vert[1];
634 sZsox = geant3->Gckine()->vert[2];
636 // Where did it hit ?
637 gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
638 gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
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;
651 // Momentum and direction of incidence
652 psx = geant3->Gctrak()->vect[6];
653 // Computing 2nd power
655 // Computing 2nd power
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();
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;
670 for (ll = 0; ll <nrooth; ++ll) rrhh[ll] = 0;
671 irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
673 ihitrak = gAlice->CurrentTrack();
675 // Flag to say that this is MIP
676 rrhh[1] = (Float_t) geant3->Gckine()->ipart;
679 rrhh[4] = (Float_t) nmodsx;
682 rrhh[6] = geant3->Gctrak()->tofg;
684 rrhh[7] = (Float_t) idpartgx;
688 // charge of current particle in electron charge unit;
689 rrhh[10] = geant3->Gckine()->charge;
691 // Zo of ckov generation
694 AddHit(ihitrak, irivol, rrhh);
699 // Current particle is cherenkov photon
700 if (geant3->Gckine()->ipart == 50) {
701 gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
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]) {
707 // fresnel losses commented out for the moment
708 // make sure that qe correction for fresnel losses is also switched off !
710 // IF (ISTOP .EQ. 2) RETURN
711 // Put transmission of electrodes in by hand
712 gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
713 cophi = TMath::Cos(TMath::ATan2(dir[0], dir[1]));
714 t = (1. - .025 / cophi) * (1. - .05 / cophi);
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;
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) {
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];
739 // copied from miphit in Freon or Quartz
740 // Where did it hit ?
741 gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
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];
747 irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
749 ihitrak = gAlice->CurrentTrack();
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];
757 // Computing 2nd power
759 // Computing 2nd power
761 rrhh[5] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
763 rrhh[6] = geant3->Gctrak()->tofg;
765 r1 = geant3->Gctrak()->upwght / 100.;
766 rrhh[7] = (Float_t) Int_t(r1+0.5);
769 rrhh[8] = geant3->Gctrak()->getot;
772 AddHit(ihitrak, irivol, rrhh);
777 if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH)
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) {
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;
789 if (geant3->Gctmed()->numed == idtmed[1003-1])
790 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=2;
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;
799 if (geant3->Gctmed()->numed == idtmed[1003-1])
800 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=12;
802 if (geant3->Gctmed()->numed == idtmed[1005-1])
803 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
805 if (geant3->Gctmed()->numed == idtmed[1009-1])
806 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
808 if (geant3->Gctmed()->numed == idtmed[1001-1])
809 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=14;
812 if (geant3->Gctmed()->numed == idtmed[1006-1])
813 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=16;
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;
827 //_____________________________________________________________________________
828 void AliRICH::RichIntegration()
831 // Integrates over RICH pads
834 TGeant3 *geant3 = (TGeant3*) gMC;
844 Int_t ixmin, ixmax, iymin, iymax;
849 Float_t y0a, qtot_check, xi1, xi2, yi1, yi2;
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
861 // IPROD = 1 PRODUCED IN FREON
862 // IPROD = 2 PRODUCED IN QUARZ
863 // get charge for MIP or cherenkov photon
865 if (geant3->Gckine()->ipart == 50) {
867 modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
871 modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
875 gMC->Gmtod(sVectIn, sVloc, 1);
876 if (TMath::Abs(sVloc[0]) >= sDlx) return;
878 if (TMath::Abs(sVloc[2]) >= sDly) return;
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;
895 for (ix = ixmin; ix <= i1; ++ix) {
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);
906 // FILL COMMON FOR PADS
913 if (sNpads > MAXSEC) {
914 // write(6,*) 'attention npads',npads
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) {
926 //isec[sNpads - 1] = sNsecon;
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
939 //isec[sNpads - 1] = sNsecon + 300;
940 } else if (ifeed == 52) {
943 //isec[sNpads - 1] = sNsecon + 200;
946 // Generate the hit for Root IO
947 for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
951 // Flag to say this is a hit
957 rrhh[6] = (Float_t) idpartgx;
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) {
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];
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];
979 // PRINT *,'RXAHIT(oldhit)=[',RRHH(20),',',RRHH(21),',',
981 // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
982 AddHit(ihitrak, irivol, rrhh);
984 for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
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;
992 rrhh[6] = (Float_t) sIqpad[sNpads - 1];
993 AddHit(ihitrak, irivol, rrhh);
998 // IF (IPART .NE. 50) WRITE(6,*) 'CHECK',QTOT,QTOT_CHECK
1000 // GENERATE FEEDBACK PHOTONS
1002 source[0] = x0 - sDlx;
1003 source[1] = sVloc[1] - .2;
1004 source[2] = y0a - sDly;
1005 gMC->Gdtom(source, source, 1);
1006 FeedBack(source, qtot);
1010 //_____________________________________________________________________________
1011 void AliRICH::AnodicWires(Float_t &y0a)
1014 // Position of anodic wires
1020 Float_t ass_i, y0, ass_i1;
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];
1036 //_____________________________________________________________________________
1037 void AliRICH::GetChargeMip(Float_t &qtot)
1040 // Calculates the charge deposited by a MIP
1051 // NUMBER OF ELECTRONS
1053 nel = (Int_t) (sDetot * 1e9 / 26.);
1057 for (i = 1; i <= i1; ++i) {
1059 qtot -= fChslope * TMath::Log(ranf[0]);
1063 //_____________________________________________________________________________
1064 void AliRICH::GetCharge(Float_t &qtot)
1074 qtot = -fChslope * TMath::Log(ranf[0]);
1077 //_____________________________________________________________________________
1078 void AliRICH::FeedBack(Float_t *source, Float_t qtot)
1081 // Generate FeedBack photons
1084 TGeant3 *geant3 = (TGeant3*) gMC;
1090 Float_t cthf, ranf[2], phif, enfp = 0, sthf;
1092 Float_t e1[3], e2[3], e3[3];
1093 Float_t vmod, uswop;
1095 Float_t dir[3], phi;
1097 Float_t pol[3], supwght;
1099 // DETERMINE NUMBER OF FEEDBACK PHOTONS
1102 fp = fAlphaFeed * qtot;
1103 nfp = gRandom->Poisson(fp);
1107 geant3->Gckin2()->ngphot = 0;
1109 for (i = 1; i <= i1; ++i) {
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);
1126 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][j]=source[j];
1129 gMC->Rndm(&random, 1);
1130 if (random <= .57) {
1132 } else if (random > .57 && random <= .7) {
1134 } else if (random > .7) {
1137 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][6] = enfp;
1138 dir[0] = sthf * TMath::Sin(phif);
1140 dir[2] = sthf * TMath::Cos(phif);
1141 gMC->Gdtom(dir, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][3], 2);
1157 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
1158 if (!vmod) for(j=0;j<3;j++) {
1164 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
1165 if (!vmod) for(j=0;j<3;j++) {
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;
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;
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;
1186 gMC->Gdtom(pol, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][7], 2);
1189 geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][10] = 0.;
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);
1196 if (geant3->Gckine()->ipart == 50 && ifeed != 52) {
1197 geant3->Gctrak()->upwght = 5100.;
1199 geant3->Gctrak()->upwght = 5200.;
1201 // geant3->Gskpho(geant3->Gckin2()->ngphot);
1202 geant3->Gctrak()->upwght = supwght;
1205 geant3->Gckin2()->ngphot = 0;
1208 //_____________________________________________________________________________
1209 Float_t AliRICH::FMathieson(Float_t lambda1, Float_t lambda2)
1212 // Mathieson charge distribution
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
1218 // K3: GEOMETRY PARAMETER
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);
1228 u1 = tanh(lambda1 * k2) * sqrtk3;
1229 u2 = tanh(lambda2 * k2) * sqrtk3;
1231 return (atan(u2) - atan(u1)) * 2 * k4;
1236 //_____________________________________________________________________________
1237 void AliRICH::UpdateMipHit(Float_t *hits)
1240 // Update some parameters when the current mip hits the detector.
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]);
1251 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
1253 //_____________________________________________________________________________
1254 void AliRICH::MakeBranch(Option_t *option)
1257 // Create a new branch in the current Root Tree.
1258 // The branch of fHits is automatically split
1260 AliDetector::MakeBranch(option);
1262 Int_t buffersize = 4000;
1263 if (gAlice->TreeH()) {
1265 gAlice->TreeH()->Branch("RICHmip",&fMips, buffersize);
1266 printf("Making Branch %s for mips\n","RICHmip");
1269 gAlice->TreeH()->Branch("RICHckov",&fCkovs, buffersize);
1270 printf("Making Branch %s for ckovs\n","RICHckov");
1273 gAlice->TreeH()->Branch("RICHpadhit",&fPadhits, buffersize);
1274 printf("Making Branch %s for padhits\n","RICHpadhit");
1279 //_____________________________________________________________________________
1280 void AliRICH::SetTreeAddress()
1283 // Set branch address for the Hits and Digits Tree.
1285 AliDetector::SetTreeAddress();
1287 TTree *treeH = gAlice->TreeH();
1290 branch = treeH->GetBranch("RICHmip");
1291 if (branch) branch->SetAddress(&fMips);
1294 branch = treeH->GetBranch("RICHckov");
1295 if (branch) branch->SetAddress(&fCkovs);
1298 branch = treeH->GetBranch("RICHpadhit");
1299 if (branch) branch->SetAddress(&fPadhits);
1304 //_____________________________________________________________________________
1305 void AliRICH::ResetHits()
1308 // Reset number of mips and the mips array for RICH
1310 AliDetector::ResetHits();
1312 if (fMips) fMips->Clear();
1313 // Reset number of Ckovs and the Ckovs array for RICH
1315 if (fCkovs) fCkovs->Clear();
1316 // Reset number of Padhits and the Padhits array for RICH
1318 if (fPadhits) fPadhits->Clear();
1321 //_____________________________________________________________________________
1322 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1325 // Distance from mouse RICH on the display
1331 //_____________________________________________________________________________
1332 void AliRICH::PreTrack()
1335 // Called before transporting a track
1346 //_____________________________________________________________________________
1347 void AliRICH::Streamer(TBuffer &R__b)
1350 // Stream an object of class AliRICH.
1352 if (R__b.IsReading()) {
1353 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1354 AliDetector::Streamer(R__b);
1358 R__b >> fMips; //diff
1359 R__b >> fCkovs; //diff
1360 R__b >> fPadhits; //diff
1366 R__b.WriteVersion(AliRICH::IsA());
1367 AliDetector::Streamer(R__b);
1371 R__b << fMips; //diff
1372 R__b << fCkovs; //diff
1373 R__b << fPadhits; //diff
1383 ///////////////////////////////////////////////////////////////////////////////
1385 // Ring Imaging Cherenkov //
1386 // This class contains version 1 of the Ring Imaging Cherenkov //
1390 <img src="picts/AliRICHv1Class.gif">
1394 ///////////////////////////////////////////////////////////////////////////////
1396 //_____________________________________________________________________________
1397 AliRICHv1::AliRICHv1() : AliRICH()
1400 // Default Constructor RICH for version 1
1404 //_____________________________________________________________________________
1405 AliRICHv1::AliRICHv1(const char *name, const char *title)
1406 : AliRICH(name,title)
1409 // Standard constructor for RICH version 1
1413 //_____________________________________________________________________________
1414 AliRICHv1::~AliRICHv1()
1417 // Standard destructor for RICH version 1
1421 //_____________________________________________________________________________
1422 void AliRICHv1::CreateGeometry()
1425 // Create the geometry for RICH version 1
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)
1433 <img src="picts/AliRICHv1.gif">
1438 <img src="picts/AliRICHv1Tree.gif">
1443 Int_t *idtmed = fIdtmed->GetArray()-999;
1450 // --- Define the RICH detector
1451 // External aluminium box
1455 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
1457 // Sensitive part of the whole RICH
1461 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
1467 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
1473 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
1479 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
1481 // Spacers (cylinders)
1485 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
1491 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
1493 // Frame of opaque quartz
1497 gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
1499 // Little bar of opaque quartz
1503 gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
1509 gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
1515 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
1521 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1527 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1533 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1535 // --- Places the detectors defined with GSVOLU
1536 // Place material inside RICH
1537 gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
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");
1544 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1546 for (i = 1; i <= 9; ++i) {
1547 zs = (5 - i) * 14.4;
1548 gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
1550 for (i = 10; i <= 18; ++i) {
1551 zs = (14 - i) * 14.4;
1552 gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
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");
1566 // Place RICH inside ALICE apparatus
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.);
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");
1586 //_____________________________________________________________________________
1587 void AliRICHv1::DrawModule()
1590 // Draw a shaded view of the Ring Imaging Cherenkov
1593 TGeant3 *geant3 = (TGeant3*) gMC;
1595 // Set everything unseen
1596 gMC->Gsatt("*", "seen", -1);
1598 // Set ALIC mother transparent
1599 gMC->Gsatt("ALIC","SEEN",0);
1601 // Set the volumes visible
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);
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();
1629 gMC->Gdraw("alic", 60, 50, 0, 10, 0, .03, .03);
1630 geant3->Gdhead(1111, "Ring Imaging Cherenkov version 1");
1631 geant3->Gdman(16, 6, "MAN");
1632 geant3->Gdopt("hide", "off");
1635 //_____________________________________________________________________________
1636 void AliRICHv1::CreateMaterials()
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)
1645 Int_t ISXFLD = gAlice->Field()->Integ();
1646 Float_t SXMGMX = gAlice->Field()->Max();
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,
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. };
1674 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1678 Int_t nlmatmet, nlmatqua;
1679 Float_t wmatquao[2], rindex_freon[14];
1681 Float_t aquao[2], epsil, stmin, zquao[2];
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];
1686 Int_t *idtmed = fIdtmed->GetArray()-999;
1688 TGeant3 *geant3 = (TGeant3*) gMC;
1690 // --- Photon energy (GeV)
1691 // --- Refraction indexes
1692 for (i = 0; i < 14; ++i) {
1693 rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
1695 // need to be changed
1697 // --- Absorbtion lenghts (in cm)
1698 // DATA ABSCO_QUARZ /
1700 // & 29.85, 7.34, 4.134, 1.273, 0.722,
1701 // & 0.365, 0.365, 0.365, 0. /
1702 // need to be changed
1704 // --- Detection efficiencies (quantum efficiency for CsI)
1705 // --- Define parameters for honeycomb.
1706 // Used carbon of equivalent rad. lenght
1713 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1724 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1735 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1746 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1757 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1764 // --- Parameters to include in GSMATE related to aluminium sheet
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);
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);
1800 // Switch on delta-ray production in the methane and freon gaps
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.);
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);
1829 ClassImp(AliRICHhit)
1831 //_____________________________________________________________________________
1832 AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol,
1833 Float_t *hits, Int_t fNpadhits) :
1837 // Standard constructor for RICH hit
1840 for (i=0;i<2;i++) fVolume[i] = vol[i];
1841 fTrack = (int) track;
1843 fPart = (int) hits[ 1];
1844 //AliHit information, position of the hit
1848 fFirstpad = fNpadhits;
1849 fLastpad = fNpadhits-1;
1852 fModule = (int) hits[ 4];
1854 fArrivaltime= hits[ 6];
1855 fFeed = (int) hits[ 7];
1858 ClassImp(AliRICHmip)
1860 //_____________________________________________________________________________
1861 AliRICHmip::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)
1866 // Standard constructor for RICH Mip hit
1873 if ((int) hits[12]){
1874 fFirstCkov = fNckovs;
1875 fLastCkov = fNckovs-1;
1882 ClassImp(AliRICHckov)
1884 //_____________________________________________________________________________
1885 AliRICHckov::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)
1890 // Standard creator for RICH Cherenkov information
1893 fStop = (int) hits[9];
1899 ClassImp(AliRICHpadhit)
1902 //_____________________________________________________________________________
1903 AliRICHpadhit::AliRICHpadhit(Int_t shunt, Int_t track, Int_t *vol,
1904 Float_t *hits, Int_t fNmips, Int_t fNckovs):
1908 // Standard constructor for RICH pad hit
1911 for (i=0;i<2;i++) fVolume[i] = vol[i];
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];
1921 fZ = -999.0; // Not implemented