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