[u/mrichter/AliRoot.git] / RICH / AliRICHv1.cxx

/**************************************************************************
 * 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 <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;
  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->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];

    
  id=gMC->CurrentVolID(copy); iCurrentChamber=copy;
  Float_t cherenkovLoss=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
    
    /********************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->GetMCApp()->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*********************
    
  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->GetMCApp()->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->GetMCApp()->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()
Commit	Line	Data
c28632f0	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
c28632f0	16
237c933d	17	#include <TParticle.h>
88cb7938	18	#include <TRandom.h>
88cb7938	19	#include <TVirtualMC.h>
d128c9d1	20	#include <TPDGCode.h>
c28632f0	21
88cb7938	22	#include "AliRICHv1.h"
c60862bf	23	#include "AliRICHParam.h"
	24	#include <AliConst.h>
	25	#include <AliPDG.h>
	26	#include <AliRun.h>
5d12ce38	27	#include <AliMC.h>
c28632f0	28
d128c9d1	29	ClassImp(AliRICHv1)
d128c9d1	30	//______________________________________________________________________________
c28632f0	31	AliRICHv1::AliRICHv1(const char name, const char title)
d128c9d1	32	:AliRICH(name,title)
3ea9cb08	33	{// Full version of RICH with hits and diagnostics
	34	if(GetDebug())Info("named ctor","Start.");
	35	if(GetDebug())Info("named ctor","Stop.");
	36	}//named ctor
	37	//______________________________________________________________________________
237c933d	38	void AliRICHv1::Init()
3ea9cb08	39	{//nothing to do here, all the work in ctor or CreateGeometry()
	40	if(GetDebug())Info("Init","Start.");
	41	if(GetDebug())Info("Init","Stop.");
	42	}//void AliRICHv1::Init()
d128c9d1	43	//______________________________________________________________________________
	44	void AliRICHv1::StepManager()
	45	{//Full Step Manager
	46
c60862bf	47	Int_t copy, id;
	48	static Int_t iCurrentChamber;
	49	static Int_t vol[2];
	50	Int_t ipart;
	51	static Float_t hits[22];
	52	static Float_t ckovData[19];
	53	TLorentzVector x4,p4;
	54	Float_t pos[3],mom[4],localPos[3],localMom[4];
	55	Float_t theta,phi,localTheta,localPhi;
	56	Float_t destep, step;
	57	Double_t ranf[2];
	58	Int_t nPads=kBad;
	59	Float_t coscerenkov;
	60	static Float_t eloss, tlength;
	61	const Float_t kBig=1.e10;
d128c9d1	62
5d12ce38	63	TParticle current = (TParticle)(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
d128c9d1	64
d128c9d1	65
d128c9d1	66
c60862bf	67	id=gMC->CurrentVolID(copy); iCurrentChamber=copy;
c60862bf	68	Float_t cherenkovLoss=0;
d128c9d1	69
c60862bf	70	gMC->TrackPosition(x4); pos[0]=x4(0); pos[1]=x4(1); pos[2]=x4(2);
	71	ckovData[1] = pos[0]; ckovData[2] = pos[1]; ckovData[3] = pos[2];
	72	ckovData[6] = 0; // dummy track length
d128c9d1	73
d128c9d1	74	/******************Store production parameters for Cerenkov photons**********************/
853634d3	75
c60862bf	76	if(gMC->TrackPid()==kCerenkov){//C
	77	Float_t ckovEnergy = current->Energy();
	78	if(ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 ){//C+E
	79	if(gMC->IsTrackEntering()){//C+E+enter
	80	if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)\|\|gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//C+E+enter+FRE
	81	if(gMC->IsNewTrack()){//C+E+enter+FRE+new
	82	Int_t mother=current->GetFirstMother();
	83	ckovData[10]=mother;
5d12ce38	84	ckovData[11]=gAlice->GetMCApp()->GetCurrentTrackNumber();
c60862bf	85	ckovData[12]=1; //Media where photon was produced 1->Freon, 2->Quarz
	86	fCkovNumber++;
	87	fFreonProd=1;
	88	}//C+E+enter+FRE+new
	89	}//C+E+enter+FRE
	90	if(gMC->IsNewTrack()&&gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) ckovData[12]=2;
	91	}//C+E+enter
	92	if(ckovData[12]==1){//C+E+produced in Freon
	93	if(gMC->IsTrackEntering()){ //is track entering?
	94	if (gMC->VolId("META")==gMC->CurrentVolID(copy)){ //is it in gap?
	95	gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
d128c9d1	96
	97	gMC->Gmtod(mom,localMom,2);
	98	Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
	99	Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
	100	/************** Photons lost in second grid have to be calculated by hand**********/
	101	gMC->GetRandom()->RndmArray(1,ranf);
	102	if (ranf[0] > t) {
	103	gMC->StopTrack();
	104	ckovData[13] = 5;
543d5224	105	//AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1	106	}
	107	/**********************************************************************************/
	108	} //gap
	109
c60862bf	110	if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){ //is it in csi?
	111
	112	gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
d128c9d1	113
	114	gMC->Gmtod(mom,localMom,2);
	115	/******* Photons lost by Fresnel reflection have to be calculated by hand******/
	116	/*********************Cerenkov phtons (always polarised)***********************/
	117	Double_t localTc = localMom[0]localMom[0]+localMom[2]localMom[2];
	118	Double_t localRt = TMath::Sqrt(localTc);
	119	localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
	120	Double_t cotheta = TMath::Abs(cos(localTheta));
	121	Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
	122	gMC->GetRandom()->RndmArray(1,ranf);
	123	if (ranf[0] < t) {
	124	gMC->StopTrack();
	125	ckovData[13] = 6;
543d5224	126	//AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1	127
	128	//printf("Added One (2)!\n");
	129	//printf("Lost by Fresnel\n");
	130	}
	131	/**********************************************************************************/
	132	}
	133	} //track entering?
d128c9d1	134	/******************Evaluation of losses**********************/
d128c9d1	135	/****************still in the old fashion********************/
d128c9d1	136	TArrayI procs;
	137	Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
	138	for (Int_t i = 0; i < i1; ++i) {
	139	// Reflection loss
	140	if (procs[i] == kPLightReflection) { //was it reflected
	141	ckovData[13]=10;
	142	if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) \|\| gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
	143	ckovData[13]=1;
	144	if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
	145	ckovData[13]=2;
d128c9d1	146	} //reflection question
	147
	148	// Absorption loss
	149	else if (procs[i] == kPLightAbsorption) { //was it absorbed?
	150	//printf("Got in absorption\n");
	151	ckovData[13]=20;
	152	if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) \|\| gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
	153	ckovData[13]=11;
	154	if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
	155	ckovData[13]=12;
	156	if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
	157	ckovData[13]=13;
	158	if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
	159	ckovData[13]=13;
	160
	161	if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
	162	ckovData[13]=15;
	163
	164	// CsI inefficiency
	165	if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
	166	ckovData[13]=16;
	167	}
	168	gMC->StopTrack();
543d5224	169	//AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1	170	} //absorption question
	171
	172
	173	// Photon goes out of tracking scope
	174	else if (procs[i] == kPStop) { //is it below energy treshold?
	175	ckovData[13]=21;
	176	gMC->StopTrack();
543d5224	177	//AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1	178	} // energy treshold question
	179	} //number of mechanisms cycle
	180	/********************End of evaluation**********************/
c60862bf	181	}//C+E+produced in Freon
	182	}//C+E
	183	}//C
	184	//***********************************End of Production Parameters Storing*******************
d128c9d1	185
c60862bf	186	if(gMC->TrackPid()==kCerenkov\|\|gMC->TrackPid()==kFeedback){//CF
	187	if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){//CF+CSI
	188	if(gMC->Edep()>0.){//CF+CSI+DE
	189	gMC->TrackPosition(x4); pos[0]=x4(0); pos[1]=x4(1); pos[2]=x4(2);
	190	gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
	191
	192	Double_t tc = mom[0]mom[0]+mom[1]mom[1];
	193	Double_t rt = TMath::Sqrt(tc);
	194	theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
	195	phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
d128c9d1	196
c60862bf	197	gMC->CurrentVolOffID(2,copy); vol[0]=copy; iCurrentChamber=vol[0];
	198
	199	gMC->Gmtod(pos,localPos,1); gMC->Gmtod(mom,localMom,2);
d128c9d1	200
c60862bf	201	Param()->SigGenInit(localPos[0],localPos[2]);
	202	ckovData[0]=gMC->TrackPid();
	203	ckovData[1]=pos[0]; ckovData[2]=pos[1]; ckovData[3]=pos[2];
	204	ckovData[4]=theta; ckovData[5]=phi; //theta-phi angles of incidence
	205	ckovData[8]=(Float_t) fNsdigits; // first sdigit
	206	ckovData[9]=-1; // last pad hit
	207	ckovData[13]=4; // photon was detected
	208	ckovData[14]=mom[0]; ckovData[15]=mom[1]; ckovData[16]=mom[2];
d128c9d1	209
c60862bf	210	destep = gMC->Edep();
	211	gMC->SetMaxStep(kBig);
	212	cherenkovLoss += destep;
	213	ckovData[7]=cherenkovLoss;
d128c9d1	214
c60862bf	215	GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
d128c9d1	216
c60862bf	217	if (fNsdigits > (Int_t)ckovData[8]) {
	218	ckovData[8]= ckovData[8]+1;
	219	ckovData[9]= (Float_t) fNsdigits;
	220	}
d128c9d1	221
c60862bf	222	ckovData[17] = nPads;
	223	AliRICHhit mipHit = (AliRICHhit) (fHits->UncheckedAt(0));
	224	if(mipHit){
	225	mom[0] = current->Px(); mom[1] = current->Py(); mom[2] = current->Pz();
	226	Float_t mipPx = mipHit->MomX(); Float_t mipPy = mipHit->MomY(); Float_t mipPz = mipHit->MomZ();
d128c9d1	227
c60862bf	228	Float_t r = mom[0]mom[0] + mom[1]mom[1] + mom[2]*mom[2];
	229	Float_t rt = TMath::Sqrt(r);
	230	Float_t mipR = mipPxmipPx + mipPymipPy + mipPz*mipPz;
	231	Float_t mipRt = TMath::Sqrt(mipR);
	232	if((rt*mipRt) > 0)
	233	coscerenkov = (mom[0]mipPx + mom[1]mipPy + mom[2]mipPz)/(rtmipRt);
	234	else
	235	coscerenkov = 0;
	236	ckovData[18]=TMath::ACos(coscerenkov);//Cerenkov angle
	237	}
5d12ce38	238	AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);//PHOTON HIT CF+CSI+DE
543d5224	239	//AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
c60862bf	240	}//CF+CSI+DE
	241	}//CF+CSI
	242	}/CF/else if(gMC->TrackCharge()){//MIP
	243	if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)\|\|gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//MIP+FRE
	244	gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
	245	hits[19]=mom[0]; hits [20] = mom[1]; hits [21] = mom[2]; fFreonProd=1;
	246	}//MIP+FRE
	247	if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)){//MIP+GAP
	248	gMC->CurrentVolOffID(3,copy); vol[0]=copy; iCurrentChamber=vol[0];
	249	gMC->TrackPosition(x4); pos[0]=x4(0); pos[1]=x4(1); pos[2]=x4(2);
	250	gMC->TrackMomentum(p4); mom[0]=p4(0); mom[1]=p4(1); mom[2]=p4(2); mom[3]=p4(3);
	251	gMC->Gmtod(pos,localPos,1); gMC->Gmtod(mom,localMom,2);
	252	ipart =gMC->TrackPid();
	253	destep = gMC->Edep();step = gMC->TrackStep();// momentum loss and steplength in last step
	254	if(gMC->IsTrackEntering()){//MIP+GAP+Enter record hit when mip enters ...
	255	Double_t tc = mom[0]mom[0]+mom[1]mom[1];
	256	Double_t rt = TMath::Sqrt(tc);
	257	theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
	258	phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
	259	Double_t localTc = localMom[0]localMom[0]+localMom[2]localMom[2];
	260	Double_t localRt = TMath::Sqrt(localTc);
	261	localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
	262	localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
	263	hits[0] = Float_t(ipart); // particle type
	264	hits[1] = localPos[0]; hits[2] = localPos[1]; hits[3] = localPos[2];
	265	hits[4] = localTheta; hits[5] = localPhi; // theta-phi angles of incidence
	266	hits[8] = (Float_t) fNsdigits; // first sdigit
	267	hits[9] = -1; // last pad hit
	268	hits[13] = fFreonProd; // did id hit the freon?
	269	hits[14] = mom[0]; hits[15] = mom[1]; hits[16] = mom[2];
	270	hits[18] = 0; // dummy cerenkov angle
	271	tlength = 0; eloss = 0; fFreonProd = 0;
	272	C(iCurrentChamber)->LocaltoGlobal(localPos,hits+1);
	273
	274	Param()->SigGenInit(localPos[0], localPos[2]);
	275	}/MIP+GAP+Enter/else if(gMC->IsTrackExiting()\|\|gMC->IsTrackStop()\|\|gMC->IsTrackDisappeared()){//MIP+GAP+Exit
	276	gMC->SetMaxStep(kBig);
	277	eloss += destep;
	278	tlength += step;
	279	if (eloss > 0) {
	280	if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
	281	GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
	282	hits[17] = nPads;
	283	}
	284	hits[6]=tlength; hits[7]=eloss;
	285	if(fNsdigits > (Int_t)hits[8]) {
	286	hits[8]= hits[8]+1;
	287	hits[9]= (Float_t) fNsdigits;
	288	}
5d12ce38	289	AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);//MIP HIT MIP+GAP+Exit
c60862bf	290	eloss = 0;
	291	}/MIP+GAP+Exit/else if(Param()->SigGenCond(localPos[0], localPos[2])){//MIP+GAP+Spec
	292	Param()->SigGenInit(localPos[0], localPos[2]);
	293	if(eloss>0){
	294	if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
	295	GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Spec
	296	hits[17] = nPads;
	297	}
	298
	299	eloss = destep; tlength += step ;
	300	}/MIP+GAP+Spec/else{//MIP+GAP+nothing special
	301	eloss += destep;
	302	tlength += step ;
	303	}//MIP+GAP+nothing special
	304	}//MIP+GAP
	305	}//MIP
d128c9d1	306	}//void AliRICHv1::StepManager()