#include <AliMC.h>
ClassImp(AliRICHv1)
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
void AliRICHv1::StepManager()
{
-//Full Step Manager
-
+// Full Step Manager.
+// 3- Ckovs absorbed on Collection electrods
+// 5- Ckovs absorbed on Cathode wires
+// 6- Ckovs absorbed on Anod wires
+
Int_t copy;
static Int_t iCurrentChamber;
- TLorentzVector x4,p4;
- Float_t pos[3],mom[4],localPos[3],localMom[4];
- Float_t coscerenkov;
-
- TParticle *current = (TParticle*)(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
-
-
-
- Float_t cherenkovLoss=0;
-
-
-
- if(gMC->TrackPid()==kCerenkov){//C
- Float_t ckovEnergy = current->Energy();
- if(ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 ){//C+E
- if(gMC->IsTrackEntering()){ //is track entering?
-
- if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){ //is it in csi?
- gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
- gMC->Gmtod(mom,localMom,2);
- Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
- Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
- Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
- if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)) gMC->StopTrack();
-
- }//C+E+produced in Freon
- } //track entering?
- }//C+E
- }//C
-
- if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI
- if(gMC->Edep()>0.){//CF+CSI+DE
- gMC->TrackPosition(x4); pos[0]=x4(0); pos[1]=x4(1); pos[2]=x4(2);
- gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
- if(IsFresnelLoss()){ gMC->StopTrack(); return;}
- gMC->CurrentVolOffID(2,copy);iCurrentChamber=copy;
+//history of Cerenkovs
+ if(gMC->TrackPid()==kCerenkov){
+ if( gMC->IsNewTrack() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(0)++;// 0- Ckovs produced in radiator
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(1)++;// 1- Ckovs absorbed in radiator
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRWI")) fCounters(2)++;// 2- Ckovs absorbed in radiator window
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RICH")) fCounters(4)++;// 4- Ckovs absorbed in CH4
+ }
+
+//Treat photons
+ static TLorentzVector cerX4;
+ if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("RPC ")){//photon in PC
+ if(gMC->Edep()>0){//photon in PC +DE
+ if(IsLostByFresnel()){
+ if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI
+ gMC->StopTrack();
+ return;
+ }
+ gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
- gMC->Gmtod(pos,localPos,1); gMC->Gmtod(mom,localMom,2);
-
- cherenkovLoss += gMC->Edep();
-
- AliRICHhit *mipHit = (AliRICHhit*) (fHits->UncheckedAt(0));
- if(mipHit){
- mom[0] = current->Px(); mom[1] = current->Py(); mom[2] = current->Pz();
- Float_t mipPx = mipHit->MomX(); Float_t mipPy = mipHit->MomY(); Float_t mipPz = mipHit->MomZ();
-
- Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
- Float_t rt = TMath::Sqrt(r);
- Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
- Float_t mipRt = TMath::Sqrt(mipR);
- if((rt*mipRt) > 0)
- coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
- else
- coscerenkov = 0;
- }
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),x4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
- GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
- }//CF+CSI+DE
- }//CF in CSI
+ AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
+ fCounters(8)++;//4- Ckovs converted to electron on CsI
+ GenerateFeedbacks(iCurrentChamber);
+ }//photon in PC and DE >0
+ }//photon in PC
//Treat charged particles
static Float_t eloss;
- if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP
- gMC->CurrentVolOffID(3,copy); iCurrentChamber=copy;
- if(gMC->IsTrackEntering()||gMC->IsNewTrack())//MIP in GAP Entering
+ static TLorentzVector mipInX4,mipOutX4;
+ if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("RGAP")){//MIP in GAP
+ gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
+ if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
eloss=0;
- else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP Exiting
+ gMC->TrackPosition(mipInX4);
+ }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared
eloss+=gMC->Edep();//take into account last step dEdX
- gMC->TrackPosition(x4);
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),x4.Vect(),eloss);//HIT for MIP for conditions: MIP in GAP Exiting
+ gMC->TrackPosition(mipOutX4);
+ AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting
GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
}else//MIP in GAP going inside
eloss += gMC->Edep();
}//MIP in GAP
-}//void AliRICHv1::StepManager()
-
-Bool_t AliRICHv1::IsFresnelLoss()
+}//StepManager()
+//__________________________________________________________________________________________________
+Bool_t AliRICHv1::IsLostByFresnel()
{
- return kFALSE;
-}
+// Calculate probability for the photon to be lost by Fresnel reflection.
+ TLorentzVector p4;
+ Double_t mom[3],localMom[3];
+ gMC->TrackMomentum(p4); mom[0]=p4(1); mom[1]=p4(2); mom[2]=p4(3);
+ localMom[0]=0; localMom[1]=0; localMom[2]=0;
+ gMC->Gmtod(mom,localMom,2);
+ Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
+ Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
+ Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
+ if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
+ AliDebug(1,"Photon lost");
+ return kTRUE;
+ }else
+ return kFALSE;
+}//IsLostByFresnel()
+//__________________________________________________________________________________________________
+void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
+{
+// Generate FeedBack photons for the current particle. To be invoked from StepManager().
+// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
+
+ TLorentzVector x4;
+ gMC->TrackPosition(x4);
+ TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
+ TVector2 xspe=x2;
+ Int_t sector=P()->Loc2Sec(xspe); if(sector==kBad) return; //hit in dead zone, nothing to produce
+ Int_t iTotQdc=P()->TotQdc(x2,eloss);
+ Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc);
+ AliDebug(1,Form("N photons=%i",iNphotons));
+ Int_t j;
+ Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
+//Generate photons
+ for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
+ Double_t ranf[2];
+ gMC->GetRandom()->RndmArray(2,ranf); //Sample direction
+ cthf=ranf[0]*2-1.0;
+ if(cthf<0) continue;
+ sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
+ phif = ranf[1] * 2 * TMath::Pi();
+
+ if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
+ enfp = 7.5e-9;
+ else if(randomNumber<=0.7)
+ enfp = 6.4e-9;
+ else
+ enfp = 7.9e-9;
+
+
+ dir[0] = sthf * TMath::Sin(phif); dir[1] = cthf; dir[2] = sthf * TMath::Cos(phif);
+ gMC->Gdtom(dir, mom, 2);
+ mom[0]*=enfp; mom[1]*=enfp; mom[2]*=enfp;
+ mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
+
+ // Polarisation
+ e1[0]= 0; e1[1]=-dir[2]; e1[2]= dir[1];
+ e2[0]=-dir[1]; e2[1]= dir[0]; e2[2]= 0;
+ e3[0]= dir[1]; e3[1]= 0; e3[2]=-dir[0];
+
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+ if (!vmod) for(j=0;j<3;j++) {
+ uswop=e1[j];
+ e1[j]=e3[j];
+ e3[j]=uswop;
+ }
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+ if (!vmod) for(j=0;j<3;j++) {
+ uswop=e2[j];
+ e2[j]=e3[j];
+ e3[j]=uswop;
+ }
+
+ 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;
+ 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;
+
+ phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
+ for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
+ gMC->Gdtom(pol, pol, 2);
+ Int_t outputNtracksStored;
+ gAlice->GetMCApp()->PushTrack(1, //do not transport
+ gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
+ kFeedback, //PID
+ mom[0],mom[1],mom[2],mom[3], //track momentum
+ x4.X(),x4.Y(),x4.Z(),x4.T(), //track origin
+ pol[0],pol[1],pol[2], //polarization
+ kPFeedBackPhoton,
+ outputNtracksStored,
+ 1.0);
+ }//feedbacks loop
+ AliDebug(1,"Stop.");
+}//GenerateFeedbacks()
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff