Using static variables to access the volume ID avoiding search by name. This reduces...
[u/mrichter/AliRoot.git] / RICH / AliRICHv1.cxx
1 // **************************************************************************
2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 // *                                                                        *
4 // * Author: The ALICE Off-line Project.                                    *
5 // * Contributors are mentioned in the code where appropriate.              *
6 // *                                                                        *
7 // * Permission to use, copy, modify and distribute this software and its   *
8 // * documentation strictly for non-commercial purposes is hereby granted   *
9 // * without fee, provided that the above copyright notice appears in all   *
10 // * copies and that both the copyright notice and this permission notice   *
11 // * appear in the supporting documentation. The authors make no claims     *
12 // * about the suitability of this software for any purpose. It is          *
13 // * provided "as is" without express or implied warranty.                  *
14 // **************************************************************************
15
16
17 #include "AliRICHv1.h"
18 #include "AliRICHParam.h"
19 #include "AliRICHChamber.h"
20 #include <TParticle.h> 
21 #include <TRandom.h> 
22 #include <TVirtualMC.h>
23 #include <TPDGCode.h>
24
25 #include <AliConst.h>
26 #include <AliPDG.h>
27 #include <AliRun.h>
28 #include <AliMC.h>
29
30 ClassImp(AliRICHv1)    
31 //__________________________________________________________________________________________________
32 void AliRICHv1::StepManager()
33 {
34 // Full Step Manager.
35 // 3- Ckovs absorbed on Collection electrods
36 // 5- Ckovs absorbed on Cathode wires
37 // 6- Ckovs absorbed on Anod wires
38          
39   Int_t          copy;
40   static Int_t   iCurrentChamber;
41   static Int_t idRRAD = gMC->VolId("RRAD");
42   static Int_t idRRWI = gMC->VolId("RRWI");
43   static Int_t idRICH = gMC->VolId("RICH");
44   static Int_t idRPC  = gMC->VolId("RPC ");
45   static Int_t idRGAP = gMC->VolId("RGAP");
46 //history of Cerenkovs
47   if(gMC->TrackPid()==kCerenkov){
48     if( gMC->IsNewTrack()   && gMC->CurrentVolID(copy)==idRRAD) fCounters(0)++;// 0- Ckovs produced in radiator
49     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRAD) fCounters(1)++;// 1- Ckovs absorbed in radiator
50     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRWI) fCounters(2)++;// 2- Ckovs absorbed in radiator window
51     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRICH) fCounters(4)++;// 4- Ckovs absorbed in CH4
52   }
53           
54 //Treat photons    
55   static TLorentzVector cerX4;
56   if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==idRPC){//photon in PC
57     if(gMC->Edep()>0){//photon in PC +DE
58       if(IsLostByFresnel()){ 
59         if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI
60         gMC->StopTrack();
61         return;
62       }        
63       gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
64         
65       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
66       fCounters(8)++;//4- Ckovs converted to electron on CsI
67       GenerateFeedbacks(iCurrentChamber);
68     }//photon in PC and DE >0 
69   }//photon in PC
70   
71 //Treat charged particles  
72   static Float_t eloss;
73   static TLorentzVector mipInX4,mipOutX4;
74   if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==idRGAP){//MIP in GAP
75     gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
76     if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
77       eloss=0;                                                           
78       gMC->TrackPosition(mipInX4);
79     }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared
80       eloss+=gMC->Edep();//take into account last step dEdX
81       gMC->TrackPosition(mipOutX4);  
82       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting
83       GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
84     }else//MIP in GAP going inside
85       eloss   += gMC->Edep();
86   }//MIP in GAP
87 }//StepManager()
88 //__________________________________________________________________________________________________
89 Bool_t AliRICHv1::IsLostByFresnel()
90 {
91 // Calculate probability for the photon to be lost by Fresnel reflection.
92   TLorentzVector p4;
93   Double_t mom[3],localMom[3];
94   gMC->TrackMomentum(p4);   mom[0]=p4(1);   mom[1]=p4(2);   mom[2]=p4(3);
95   localMom[0]=0; localMom[1]=0; localMom[2]=0;
96   gMC->Gmtod(mom,localMom,2);
97   Double_t localTc    = localMom[0]*localMom[0]+localMom[2]*localMom[2];
98   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
99   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
100   if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
101     AliDebug(1,"Photon lost");
102     return kTRUE;
103   }else
104     return kFALSE;
105 }//IsLostByFresnel()
106 //__________________________________________________________________________________________________
107 void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
108 {
109 // Generate FeedBack photons for the current particle. To be invoked from StepManager().
110 // eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
111   
112   TLorentzVector x4;
113   gMC->TrackPosition(x4);  
114   TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
115   TVector2 xspe=x2;
116   Int_t sector=P()->Loc2Sec(xspe);  if(sector==kBad) return; //hit in dead zone, nothing to produce
117   Int_t iTotQdc=P()->TotQdc(x2,eloss);
118   Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc);    
119   AliDebug(1,Form("N photons=%i",iNphotons));
120   Int_t j;
121   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
122 //Generate photons
123   for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
124     Double_t ranf[2];
125     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
126     cthf=ranf[0]*2-1.0;
127     if(cthf<0) continue;
128     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
129     phif = ranf[1] * 2 * TMath::Pi();
130     
131     if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
132       enfp = 7.5e-9;
133     else if(randomNumber<=0.7)
134       enfp = 6.4e-9;
135     else
136       enfp = 7.9e-9;
137     
138
139     dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
140     gMC->Gdtom(dir, mom, 2);
141     mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
142     mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
143     
144     // Polarisation
145     e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
146     e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
147     e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
148     
149     vmod=0;
150     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
151     if (!vmod) for(j=0;j<3;j++) {
152       uswop=e1[j];
153       e1[j]=e3[j];
154       e3[j]=uswop;
155     }
156     vmod=0;
157     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
158     if (!vmod) for(j=0;j<3;j++) {
159       uswop=e2[j];
160       e2[j]=e3[j];
161       e3[j]=uswop;
162     }
163     
164     vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;    
165     vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
166     
167     phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
168     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
169     gMC->Gdtom(pol, pol, 2);
170     Int_t outputNtracksStored;    
171     gAlice->GetMCApp()->PushTrack(1,                 //do not transport
172                      gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track 
173                      kFeedback,                      //PID
174                      mom[0],mom[1],mom[2],mom[3],    //track momentum  
175                      x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin 
176                      pol[0],pol[1],pol[2],           //polarization
177                      kPFeedBackPhoton,
178                      outputNtracksStored,
179                      1.0);    
180   }//feedbacks loop
181   AliDebug(1,"Stop.");
182 }//GenerateFeedbacks()