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 in Collection electrods
  36 // 5- Ckovs absorbed in Cathode wires
  37 // 6- Ckovs absorbed in Anod wires
  38
  39   Int_t          copy;
  40   static Int_t   iCurrentChamber;
  41 //history of Cerenkovs
  42   if(gMC->TrackPid()==kCerenkov){
  43     if(gMC->IsNewTrack() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(0)++;// 0- Ckovs produced in Freon
  44     if(!gMC->IsTrackAlive() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(1)++;// 1- Ckovs absorbed in Freon
  45     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("QUAR")) fCounters(2)++;//2- Ckovs absorbed in Quartz
  46     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("META")) fCounters(4)++;//4- Ckovs absorbed in CH4
  47   }
  48
  49 //Treat photons
  50   static TLorentzVector cerX4;
  51   if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI
  52     if(gMC->Edep()>0){//CF+CSI+DE
  53       if(IsLostByFresnel()){
  54         if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI
  55         gMC->StopTrack();
  56         return;
  57       }
  58       gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);
  59
  60       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
  61       fCounters(8)++;//4- Ckovs converted in CsI
  62       GenerateFeedbacks(iCurrentChamber);
  63     }//CF+CSI+DE
  64   }//CF in CSI
  65
  66 //Treat charged particles
  67   static Float_t eloss;
  68   static TLorentzVector mipInX4,mipOutX4;
  69   if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP
  70     gMC->CurrentVolOffID(3,iCurrentChamber);
  71     if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
  72       eloss=0;
  73       gMC->TrackPosition(mipInX4);
  74     }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared
  75       eloss+=gMC->Edep();//take into account last step dEdX
  76       gMC->TrackPosition(mipOutX4);
  77       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting
  78       GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
  79     }else//MIP in GAP going inside
  80       eloss   += gMC->Edep();
  81   }//MIP in GAP
  82 }//StepManager()
  83 //__________________________________________________________________________________________________
  84 Bool_t AliRICHv1::IsLostByFresnel()
  85 {
  86 // Calculate probability for the photon to be lost by Fresnel reflection.
  87   TLorentzVector p4;
  88   Double_t mom[3],localMom[3];
  89   gMC->TrackMomentum(p4);   mom[0]=p4(1);   mom[1]=p4(2);   mom[2]=p4(3);
  90   localMom[0]=0; localMom[1]=0; localMom[2]=0;
  91   gMC->Gmtod(mom,localMom,2);
  92   Double_t localTc    = localMom[0]*localMom[0]+localMom[2]*localMom[2];
  93   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
  94   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
  95   if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
  96     if(GetDebug()) Info("IsLostByFresnel","Photon lost");
  97     return kTRUE;
  98   }else
  99     return kFALSE;
 100 }//IsLostByFresnel()
 101 //__________________________________________________________________________________________________
 102 void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
 103 {
 104 // Generate FeedBack photons for the current particle. To be invoked from StepManager().
 105 // eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
 106
 107   TLorentzVector x4;
 108   gMC->TrackPosition(x4);
 109   TVector2 x2=C(iChamber)->Glob2Loc(x4);
 110   Int_t sector=P()->Loc2Sec(x2);  if(sector==kBad) return; //hit in dead zone, nothing to produce
 111   Int_t iTotQdc=P()->TotQdc(x2,eloss);
 112   Int_t iNphotons=gMC->GetRandom()->Poisson(P()->AlphaFeedback(sector)*iTotQdc);
 113   if(GetDebug())Info("GenerateFeedbacks","N photons=%i",iNphotons);
 114   Int_t j;
 115   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
 116 //Generate photons
 117   for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
 118     Double_t ranf[2];
 119     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
 120     cthf=ranf[0]*2-1.0;
 121     if(cthf<0) continue;
 122     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 123     phif = ranf[1] * 2 * TMath::Pi();
 124
 125     if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
 126       enfp = 7.5e-9;
 127     else if(randomNumber<=0.7)
 128       enfp = 6.4e-9;
 129     else
 130       enfp = 7.9e-9;
 131
 132
 133     dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
 134     gMC->Gdtom(dir, mom, 2);
 135     mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
 136     mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
 137
 138     // Polarisation
 139     e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
 140     e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
 141     e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
 142
 143     vmod=0;
 144     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 145     if (!vmod) for(j=0;j<3;j++) {
 146       uswop=e1[j];
 147       e1[j]=e3[j];
 148       e3[j]=uswop;
 149     }
 150     vmod=0;
 151     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 152     if (!vmod) for(j=0;j<3;j++) {
 153       uswop=e2[j];
 154       e2[j]=e3[j];
 155       e3[j]=uswop;
 156     }
 157
 158     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;
 159     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;
 160
 161     phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
 162     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 163     gMC->Gdtom(pol, pol, 2);
 164     Int_t outputNtracksStored;
 165     gAlice->GetMCApp()->PushTrack(1,                 //do not transport
 166                      gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
 167                      kFeedback,                      //PID
 168                      mom[0],mom[1],mom[2],mom[3],    //track momentum
 169                      x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin
 170                      pol[0],pol[1],pol[2],           //polarization
 171                      kPFeedBackPhoton,
 172                      outputNtracksStored,
 173                      1.0);
 174   }//feedbacks loop
 175   if(GetDebug()) Info("GenerateFeedbacks","Stop.");
 176 }//GenerateFeedbacks()