No return in void function (Alpha)
[u/mrichter/AliRoot.git] / RICH / AliRICHv1.cxx
index 1021965..b4f3765 100644 (file)
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
+// **************************************************************************
+// * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+// *                                                                        *
+// * Author: The ALICE Off-line Project.                                    *
+// * Contributors are mentioned in the code where appropriate.              *
+// *                                                                        *
+// * Permission to use, copy, modify and distribute this software and its   *
+// * documentation strictly for non-commercial purposes is hereby granted   *
+// * without fee, provided that the above copyright notice appears in all   *
+// * copies and that both the copyright notice and this permission notice   *
+// * appear in the supporting documentation. The authors make no claims     *
+// * about the suitability of this software for any purpose. It is          *
+// * provided "as is" without express or implied warranty.                  *
+// **************************************************************************
 
 
+#include "AliRICHv1.h"
+#include "AliRICHParam.h"
+#include "AliRICHChamber.h"
 #include <TParticle.h> 
 #include <TRandom.h> 
 #include <TVirtualMC.h>
 #include <TPDGCode.h>
 
-#include "AliRICHv1.h"
-#include "AliRICHParam.h"
 #include <AliConst.h>
 #include <AliPDG.h>
 #include <AliRun.h>
+#include <AliMC.h>
 
 ClassImp(AliRICHv1)    
-//______________________________________________________________________________
-AliRICHv1::AliRICHv1(const char *name, const char *title)
-          :AliRICH(name,title)
-{// Full version of RICH with hits and diagnostics
-  if(GetDebug())Info("named ctor","Start.");
-  if(GetDebug())Info("named ctor","Stop.");  
-}//named ctor
-//______________________________________________________________________________
-void AliRICHv1::Init()
-{//nothing to do here, all the work in ctor or CreateGeometry()
-  if(GetDebug())Info("Init","Start.");
-  if(GetDebug())Info("Init","Stop.");    
-}//void AliRICHv1::Init()
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICHv1::StepManager()
-{//Full Step Manager
-
-  Int_t          copy, id;
+{
+// 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;
-  static Int_t   vol[2];
-  Int_t          ipart;
-  static Float_t hits[22];
-  static Float_t ckovData[19];
-  TLorentzVector x4,p4;
-  Float_t        pos[3],mom[4],localPos[3],localMom[4];
-  Float_t        theta,phi,localTheta,localPhi;
-  Float_t        destep, step;
-  Double_t        ranf[2];
-  Int_t          nPads=kBad;
-  Float_t        coscerenkov;
-  static Float_t eloss,  tlength;
-  const  Float_t kBig=1.e10;
-       
-  TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->GetCurrentTrackNumber()];
+  static Int_t idRRAD = gMC->VolId("RRAD");
+  static Int_t idRRWI = gMC->VolId("RRWI");
+  static Int_t idRICH = gMC->VolId("RICH");
+  static Int_t idRPC  = gMC->VolId("RPC ");
+  static Int_t idRGAP = gMC->VolId("RGAP");
+//history of Cerenkovs
+  if(gMC->TrackPid()==kCerenkov){
+    if( gMC->IsNewTrack()   && gMC->CurrentVolID(copy)==idRRAD) fCounters(0)++;// 0- Ckovs produced in radiator
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRAD) fCounters(1)++;// 1- Ckovs absorbed in radiator
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRWI) fCounters(2)++;// 2- Ckovs absorbed in radiator window
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRICH) fCounters(4)++;// 4- Ckovs absorbed in CH4
+  }
+          
+//Treat photons    
+  static TLorentzVector cerX4;
+  if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==idRPC){//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
+       
+      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;
+  static TLorentzVector mipInX4,mipOutX4;
+  if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==idRGAP){//MIP in GAP
+    gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
+    if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
+      eloss=0;                                                           
+      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(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
+}//StepManager()
+//__________________________________________________________________________________________________
+Bool_t AliRICHv1::IsLostByFresnel()
+{
+// 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]);
     
-  id=gMC->CurrentVolID(copy); iCurrentChamber=copy;
-  Float_t cherenkovLoss=0;
+    // 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];
     
-  gMC->TrackPosition(x4);    pos[0]=x4(0);    pos[1]=x4(1);    pos[2]=x4(2);
-  ckovData[1] = pos[0];   ckovData[2] = pos[1];   ckovData[3] = pos[2]; 
-  ckovData[6] = 0;                      // dummy track length
+    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;
+    }
     
-    /********************Store production parameters for Cerenkov photons************************/ 
-
-  if(gMC->TrackPid()==kCerenkov){//C
-    Float_t ckovEnergy = current->Energy();
-    if(ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 ){//C+E
-      if(gMC->IsTrackEntering()){//C+E+enter
-        if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//C+E+enter+FRE
-          if(gMC->IsNewTrack()){//C+E+enter+FRE+new
-            Int_t mother=current->GetFirstMother(); 
-            ckovData[10]=mother;
-            ckovData[11]=gAlice->GetCurrentTrackNumber();
-            ckovData[12]=1;             //Media where photon was produced 1->Freon, 2->Quarz
-            fCkovNumber++;
-            fFreonProd=1;
-          }//C+E+enter+FRE+new
-        }//C+E+enter+FRE
-        if(gMC->IsNewTrack()&&gMC->VolId("QUAR")==gMC->CurrentVolID(copy))  ckovData[12]=2;            
-      }//C+E+enter           
-      if(ckovData[12]==1){//C+E+produced in Freon
-        if(gMC->IsTrackEntering()){                                     //is track entering?
-                   if (gMC->VolId("META")==gMC->CurrentVolID(copy)){                //is it in gap?      
-                  gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-                       
-                       gMC->Gmtod(mom,localMom,2);
-                       Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
-                       Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
-                       /**************** Photons lost in second grid have to be calculated by hand************/ 
-                       gMC->GetRandom()->RndmArray(1,ranf);
-                       if (ranf[0] > t) {
-                         gMC->StopTrack();
-                         ckovData[13] = 5;
-                         AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                       }
-                       /**********************************************************************************/
-                     }    //gap
-                   
-                   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);
-                       /********* Photons lost by Fresnel reflection have to be calculated by hand********/ 
-                       /***********************Cerenkov phtons (always polarised)*************************/
-                       Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
-                       Double_t localRt = TMath::Sqrt(localTc);
-                       localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
-                       Double_t cotheta = TMath::Abs(cos(localTheta));
-                       Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
-                           gMC->GetRandom()->RndmArray(1,ranf);
-                           if (ranf[0] < t) {
-                             gMC->StopTrack();
-                             ckovData[13] = 6;
-                             AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                               
-                             //printf("Added One (2)!\n");
-                             //printf("Lost by Fresnel\n");
-                           }
-                           /**********************************************************************************/
-                     }
-                 } //track entering?
-                 /********************Evaluation of losses************************/
-                 /******************still in the old fashion**********************/
-                 TArrayI procs;
-                 Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
-                 for (Int_t i = 0; i < i1; ++i) {
-                   //        Reflection loss 
-                   if (procs[i] == kPLightReflection) {        //was it reflected
-                     ckovData[13]=10;
-                     if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
-                       ckovData[13]=1;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
-                       ckovData[13]=2;
-                   } //reflection question
-                    
-                   //        Absorption loss 
-                   else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
-                     //printf("Got in absorption\n");
-                     ckovData[13]=20;
-                     if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
-                       ckovData[13]=11;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
-                       ckovData[13]=12;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("META")) 
-                       ckovData[13]=13;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("GAP ")) 
-                       ckovData[13]=13;
-                     
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC")) 
-                       ckovData[13]=15;
-                     
-                     //        CsI inefficiency 
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
-                       ckovData[13]=16;
-                     }
-                     gMC->StopTrack();
-                     AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                   } //absorption question 
-                   
-                   
-                   //        Photon goes out of tracking scope 
-                   else if (procs[i] == kPStop) {                 //is it below energy treshold?
-                     ckovData[13]=21;
-                     gMC->StopTrack();
-                     AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                   }   // energy treshold question         
-                 }  //number of mechanisms cycle
-                 /**********************End of evaluation************************/
-      }//C+E+produced in Freon
-    }//C+E
-  }//C
-//*************************************End of Production Parameters Storing*********************
+    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;
     
-  if(gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback){//CF      
-    if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){//CF+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);
-        
-       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
-       Double_t rt = TMath::Sqrt(tc);
-       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
-       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
-               
-       gMC->CurrentVolOffID(2,copy);   vol[0]=copy;    iCurrentChamber=vol[0];
-       
-        gMC->Gmtod(pos,localPos,1);     gMC->Gmtod(mom,localMom,2);
-
-       Param()->SigGenInit(localPos[0],localPos[2]);
-        ckovData[0]=gMC->TrackPid();   
-        ckovData[1]=pos[0];      ckovData[2]=pos[1];    ckovData[3]=pos[2];            
-        ckovData[4]=theta;       ckovData[5]=phi;       //theta-phi angles of incidence 
-        ckovData[8]=(Float_t) fNsdigits;      // first sdigit
-        ckovData[9]=-1;                       // last pad hit
-        ckovData[13]=4;                       // photon was detected
-        ckovData[14]=mom[0];     ckovData[15]=mom[1];    ckovData[16]=mom[2];
-                   
-       destep = gMC->Edep();
-       gMC->SetMaxStep(kBig);
-       cherenkovLoss  += destep;
-       ckovData[7]=cherenkovLoss;
-                   
-       GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
-                                   
-       if (fNsdigits > (Int_t)ckovData[8]) {
-         ckovData[8]= ckovData[8]+1;
-          ckovData[9]= (Float_t) fNsdigits;
-        }
-
-        ckovData[17] = nPads;
-        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;
-          ckovData[18]=TMath::ACos(coscerenkov);//Cerenkov angle
-        }
-        AddHit(gAlice->GetCurrentTrackNumber(),vol,ckovData);//PHOTON HIT CF+CSI+DE
-        AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-      }//CF+CSI+DE
-    }//CF+CSI
-  }/*CF*/else if(gMC->TrackCharge()){//MIP
-    if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//MIP+FRE
-      gMC->TrackMomentum(p4);      mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-      hits[19]=mom[0];     hits [20] = mom[1];     hits [21] = mom[2];     fFreonProd=1;
-    }//MIP+FRE
-    if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)){//MIP+GAP
-      gMC->CurrentVolOffID(3,copy);    vol[0]=copy;    iCurrentChamber=vol[0];
-      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);
-      gMC->Gmtod(pos,localPos,1);          gMC->Gmtod(mom,localMom,2);
-      ipart =gMC->TrackPid();
-      destep = gMC->Edep();step   = gMC->TrackStep();// momentum loss and steplength in last step
-      if(gMC->IsTrackEntering()){//MIP+GAP+Enter  record hit when mip enters ...
-        Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
-        Double_t rt = TMath::Sqrt(tc);
-        theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
-        phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
-       Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
-       Double_t localRt = TMath::Sqrt(localTc);
-       localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;                       
-       localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;    
-        hits[0] = Float_t(ipart);         // particle type
-        hits[1] = localPos[0];    hits[2] = localPos[1];     hits[3] = localPos[2];                 
-        hits[4] = localTheta;     hits[5] = localPhi;              // theta-phi angles of incidence 
-        hits[8] = (Float_t) fNsdigits;    // first sdigit
-        hits[9] = -1;                     // last pad hit
-        hits[13] = fFreonProd;           // did id hit the freon?
-        hits[14] = mom[0];        hits[15] = mom[1];        hits[16] = mom[2];
-        hits[18] = 0;               // dummy cerenkov angle
-        tlength = 0;        eloss   = 0;        fFreonProd = 0;
-        C(iCurrentChamber)->LocaltoGlobal(localPos,hits+1);
-          
-        Param()->SigGenInit(localPos[0], localPos[2]);
-      }/*MIP+GAP+Enter*/else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP+GAP+Exit
-        gMC->SetMaxStep(kBig);
-       eloss   += destep;
-       tlength += step;
-        if (eloss > 0) {
-          if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
-          GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
-          hits[17] = nPads;
-        }
-        hits[6]=tlength;        hits[7]=eloss;
-        if(fNsdigits > (Int_t)hits[8]) {
-          hits[8]= hits[8]+1;
-          hits[9]= (Float_t) fNsdigits;
-        }
-       AddHit(gAlice->GetCurrentTrackNumber(),vol,hits);//MIP HIT MIP+GAP+Exit
-       eloss = 0; 
-      }/*MIP+GAP+Exit*/else if(Param()->SigGenCond(localPos[0], localPos[2])){//MIP+GAP+Spec
-        Param()->SigGenInit(localPos[0], localPos[2]);
-        if(eloss>0){
-          if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
-          GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Spec
-          hits[17] = nPads;
-        }
-        
-       eloss    = destep;      tlength += step ;               
-      }/*MIP+GAP+Spec*/else{//MIP+GAP+nothing special
-        eloss   += destep;
-        tlength += step ;
-      }//MIP+GAP+nothing special
-    }//MIP+GAP
-  }//MIP
-}//void AliRICHv1::StepManager()
+    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()