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