Code cleaning, all wranings removed with new Makefile options
[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 /* $Id$ */
17
18 #include "Riostream.h"
19
20 #include <TNode.h> 
21 #include <TParticle.h> 
22 #include <TRandom.h> 
23 #include <TTUBE.h>
24 #include <TVirtualMC.h>
25 #include <TPDGCode.h>
26
27 #include "AliConst.h" 
28 #include "AliPDG.h" 
29 #include "AliRICHGeometry.h"
30 #include "AliRICHResponse.h"
31 #include "AliRICHResponseV0.h"
32 #include "AliRICHSegmentationV1.h"
33 #include "AliRICHv1.h"
34 #include "AliRun.h"
35
36 ClassImp(AliRICHv1)    
37 //______________________________________________________________________________
38 AliRICHv1::AliRICHv1(const char *name, const char *title)
39           :AliRICH(name,title)
40 {
41
42 // Full version of RICH with hits and diagnostics
43
44   // Version 0
45 // Default Segmentation, no hits
46     AliRICHSegmentationV1* segmentation = new AliRICHSegmentationV1;
47 //
48 //  Segmentation parameters
49     segmentation->SetPadSize(0.84,0.80);
50     segmentation->SetDAnod(0.84/2);
51 //
52 //  Geometry parameters
53     AliRICHGeometry* geometry = new AliRICHGeometry;
54     geometry->SetGapThickness(8);
55     geometry->SetProximityGapThickness(.4);
56     geometry->SetQuartzLength(133);
57     geometry->SetQuartzWidth(127.9);
58     geometry->SetQuartzThickness(.5);
59     geometry->SetOuterFreonLength(133);
60     geometry->SetOuterFreonWidth(41.3);
61     geometry->SetInnerFreonLength(133);
62     geometry->SetInnerFreonWidth(41.3);
63     geometry->SetFreonThickness(1.5);
64 //
65 //  Response parameters
66     AliRICHResponseV0*  response   = new AliRICHResponseV0;
67     response->SetSigmaIntegration(5.);
68     response->SetChargeSlope(27.);
69     response->SetChargeSpread(0.18, 0.18);
70     response->SetMaxAdc(4096);
71     response->SetAlphaFeedback(0.036);
72     response->SetEIonisation(26.e-9);
73     response->SetSqrtKx3(0.77459667);
74     response->SetKx2(0.962);
75     response->SetKx4(0.379);
76     response->SetSqrtKy3(0.77459667);
77     response->SetKy2(0.962);
78     response->SetKy4(0.379);
79     response->SetPitch(0.25);
80     response->SetWireSag(1);                     // 1->On, 0->Off
81     response->SetVoltage(2150);                  // Should only be 2000, 2050, 2100 or 2150
82
83 //
84 //    AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH"); 
85     
86     fCkovNumber=0;
87     fFreonProd=0;
88     Int_t i=0;
89     
90     fChambers = new TObjArray(kNCH);
91     for (i=0; i<kNCH; i++) {
92       
93       //PH      (*fChambers)[i] = new AliRICHChamber();  
94       fChambers->AddAt(new AliRICHChamber(), i);  
95       
96     }
97   
98     for (i=0; i<kNCH; i++) {
99       SetGeometryModel(i,geometry);
100       SetSegmentationModel(i, segmentation);
101       SetResponseModel(i, response);
102     }
103
104
105 }
106
107 void AliRICHv1::Init()
108 {
109
110   if(fDebug) {
111     printf("%s: *********************************** RICH_INIT ***********************************\n",ClassName());
112     printf("%s: *                                                                               *\n",ClassName());
113     printf("%s: *                        AliRICHv1 Full version started                         *\n",ClassName());
114     printf("%s: *                                                                               *\n",ClassName());
115   }
116
117   
118   AliSegmentation*  segmentation;
119   AliRICHGeometry*  geometry;
120   AliRICHResponse*  response;
121
122
123     // 
124     // Initialize Tracking Chambers
125     //
126     for (Int_t i=0; i<kNCH; i++) {
127         //printf ("i:%d",i);
128       //PH      ( (AliRICHChamber*) (*fChambers)[i])->Init(i);  
129         ( (AliRICHChamber*)fChambers->At(i))->Init(i);  
130     }  
131     
132     //
133     // Set the chamber (sensitive region) GEANT identifier
134     
135     //PH    ((AliRICHChamber*)(*fChambers)[0])->SetGid(1);  
136     //PH    ((AliRICHChamber*)(*fChambers)[1])->SetGid(2);  
137     //PH    ((AliRICHChamber*)(*fChambers)[2])->SetGid(3);  
138     //PH    ((AliRICHChamber*)(*fChambers)[3])->SetGid(4);  
139     //PH    ((AliRICHChamber*)(*fChambers)[4])->SetGid(5);  
140     //PH    ((AliRICHChamber*)(*fChambers)[5])->SetGid(6);  
141     //PH    ((AliRICHChamber*)(*fChambers)[6])->SetGid(7); 
142
143     ((AliRICHChamber*)fChambers->At(0))->SetGid(1);  
144     ((AliRICHChamber*)fChambers->At(1))->SetGid(2);  
145     ((AliRICHChamber*)fChambers->At(2))->SetGid(3);  
146     ((AliRICHChamber*)fChambers->At(3))->SetGid(4);  
147     ((AliRICHChamber*)fChambers->At(4))->SetGid(5);  
148     ((AliRICHChamber*)fChambers->At(5))->SetGid(6);  
149     ((AliRICHChamber*)fChambers->At(6))->SetGid(7);  
150
151
152     segmentation=Chamber(0).GetSegmentationModel();
153     geometry=Chamber(0).GetGeometryModel();
154     response=Chamber(0).GetResponseModel();
155     
156     Float_t offset       = 490 + 1.276 - geometry->GetGapThickness()/2;        //distance from center of mother volume to methane
157     Float_t deltaphi     = 19.5;                                               //phi angle between center of chambers - z direction
158     Float_t deltatheta   = 20;                                                 //theta angle between center of chambers - x direction
159     Float_t cosphi       = TMath::Cos(deltaphi*TMath::Pi()/180);
160     Float_t sinphi       = TMath::Sin(deltaphi*TMath::Pi()/180);
161     Float_t costheta     = TMath::Cos(deltatheta*TMath::Pi()/180);
162     Float_t sintheta     = TMath::Sin(deltatheta*TMath::Pi()/180);
163
164     Float_t pos1[3]={0.                , offset*cosphi         , offset*sinphi};
165     Float_t pos2[3]={offset*sintheta   , offset*costheta       , 0. };
166     Float_t pos3[3]={0.                , offset                , 0.};
167     Float_t pos4[3]={-offset*sintheta  , offset*costheta       , 0.};
168     Float_t pos5[3]={offset*sinphi     , offset*costheta*cosphi, -offset*sinphi};
169     Float_t pos6[3]={0.                , offset*cosphi         , -offset*sinphi};
170     Float_t pos7[3]={ -offset*sinphi   , offset*costheta*cosphi, -offset*sinphi};
171
172     Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90., 0.               , 90. - deltaphi, 90.             , deltaphi, -90.           ));
173     Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90., -deltatheta      , 90.           , 90.- deltatheta , 0.      , 0.             ));
174     Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90., 0.               , 90.           , 90.             , 0.      , 0.             ));
175     Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],new TRotMatrix("rot996","rot996",90.,  deltatheta      , 90.           , 90 + deltatheta , 0.      , 0.             ));
176     Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2         , 90.- deltatheta ,18.2     , 90 - deltatheta));
177     Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90., 0.               , 90 + deltaphi , 90.             , deltaphi, 90.            ));
178     Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90., deltatheta       , 108.2         , 90.+ deltatheta ,18.2     , 90 + deltatheta));   
179      
180     if(GetDebug()) {
181       printf("%s: *                            Pads            : %3dx%3d                          *\n",
182              ClassName(),segmentation->Npx(),segmentation->Npy());
183       printf("%s: *                            Pad size        : %5.2f x%5.2f mm2                 *\n",
184              ClassName(),segmentation->Dpx(),segmentation->Dpy()); 
185       printf("%s: *                            Gap Thickness   : %5.1f cm                         *\n",
186              ClassName(),geometry->GetGapThickness());
187       printf("%s: *                            Radiator Width  : %5.1f cm                         *\n",
188              ClassName(),geometry->GetQuartzWidth());
189       printf("%s: *                            Radiator Length : %5.1f cm                         *\n",
190              ClassName(),geometry->GetQuartzLength());
191       printf("%s: *                            Freon Thickness : %5.1f cm                         *\n",
192              ClassName(),geometry->GetFreonThickness());
193       printf("%s: *                            Charge Slope    : %5.1f ADC                        *\n",
194              ClassName(),response->ChargeSlope());
195       printf("%s: *                            Feedback Prob.  : %5.2f %%                          *\n",
196              ClassName(),response->AlphaFeedback()*100);
197       printf("%s: *                                                                               *\n",
198              ClassName());
199       printf("%s: *********************************************************************************\n",
200              ClassName());
201     }
202 }
203
204 //______________________________________________________________________________
205 void AliRICHv1::StepManager()
206 {//Full Step Manager
207
208     Int_t          copy, id;
209     static Int_t   idvol;
210     static Int_t   vol[2];
211     Int_t          ipart;
212     static Float_t hits[22];
213     static Float_t ckovData[19];
214     TLorentzVector position;
215     TLorentzVector momentum;
216     Float_t        pos[3];
217     Float_t        mom[4];
218     Float_t        localPos[3];
219     Float_t        localMom[4];
220     Float_t        localTheta,localPhi;
221     Float_t        theta,phi;
222     Float_t        destep, step;
223     Double_t        ranf[2];
224     Int_t          nPads=-1;
225     Float_t        coscerenkov;
226     static Float_t eloss, xhit, yhit, tlength;
227     const  Float_t kBig=1.e10;
228        
229     TClonesArray &lhits = *fHits;
230     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->GetCurrentTrackNumber()];
231
232     
233  
234     id=gMC->CurrentVolID(copy);
235     idvol = copy-1;
236     Float_t cherenkovLoss=0;
237     
238     gMC->TrackPosition(position);
239     pos[0]=position(0);
240     pos[1]=position(1);
241     pos[2]=position(2);
242     ckovData[1] = pos[0];                 // X-position for hit
243     ckovData[2] = pos[1];                 // Y-position for hit
244     ckovData[3] = pos[2];                 // Z-position for hit
245     ckovData[6] = 0;                      // dummy track length
246     
247     /********************Store production parameters for Cerenkov photons************************/ 
248
249     if (gMC->TrackPid() == 50000050) { //is it a Cerenkov photon? 
250
251       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
252         //{                    
253           Float_t ckovEnergy = current->Energy();
254           //energy interval for tracking
255           if  (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )       
256             //if (ckovEnergy > 0)
257             {
258               if (gMC->IsTrackEntering()){        //is track entering?
259                 //printf("Track entered (1)\n");
260                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
261                   {                                                          //is it in freo?
262                     if (gMC->IsNewTrack()){                          //is it the first step?
263                       //printf("I'm in!\n");
264                       Int_t mother = current->GetFirstMother(); 
265                       
266                       
267                       ckovData[10] = mother;
268                       ckovData[11] = gAlice->GetCurrentTrackNumber();
269                       ckovData[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
270                       //printf("Produced in FREO\n");
271                       fCkovNumber++;
272                       fFreonProd=1;
273                     }    //first step question
274                   }        //freo question
275                 
276                 if (gMC->IsNewTrack()){                                  //is it first step?
277                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
278                     {
279                       ckovData[12] = 2;
280                     }    //quarz question
281                 }        //first step question
282                 
283               }   //track entering question
284               
285               if (ckovData[12] == 1)                                        //was it produced in Freon?
286                 //if (fFreonProd == 1)
287                 {
288                   if (gMC->IsTrackEntering()){                                     //is track entering?
289                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?      
290                       {
291                         //printf("Got in META\n");
292                         gMC->TrackMomentum(momentum);
293                         mom[0]=momentum(0);
294                         mom[1]=momentum(1);
295                         mom[2]=momentum(2);
296                         mom[3]=momentum(3);
297                         
298                         gMC->Gmtod(mom,localMom,2);
299                         Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
300                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
301                         /**************** Photons lost in second grid have to be calculated by hand************/ 
302                         gMC->GetRandom()->RndmArray(1,ranf);
303                         if (ranf[0] > t) {
304                           gMC->StopTrack();
305                           ckovData[13] = 5;
306                           AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
307                         }
308                         /**********************************************************************************/
309                       }    //gap
310                     
311                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?      
312                       {
313                         //printf("Got in CSI\n");
314                         gMC->TrackMomentum(momentum);
315                         mom[0]=momentum(0);
316                         mom[1]=momentum(1);
317                         mom[2]=momentum(2);
318                         mom[3]=momentum(3);
319
320                         gMC->Gmtod(mom,localMom,2);
321                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/ 
322                         /***********************Cerenkov phtons (always polarised)*************************/
323                         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
324                         Double_t localRt = TMath::Sqrt(localTc);
325                         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
326                         Double_t cotheta = TMath::Abs(cos(localTheta));
327                         Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
328                             gMC->GetRandom()->RndmArray(1,ranf);
329                             if (ranf[0] < t) {
330                               gMC->StopTrack();
331                               ckovData[13] = 6;
332                               AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
333                                 
334                               //printf("Added One (2)!\n");
335                               //printf("Lost by Fresnel\n");
336                             }
337                             /**********************************************************************************/
338                       }
339                   } //track entering?
340                   
341                   
342                   /********************Evaluation of losses************************/
343                   /******************still in the old fashion**********************/
344                   
345                   TArrayI procs;
346                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
347                   for (Int_t i = 0; i < i1; ++i) {
348                     //        Reflection loss 
349                     if (procs[i] == kPLightReflection) {        //was it reflected
350                       ckovData[13]=10;
351                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
352                         ckovData[13]=1;
353                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
354                         ckovData[13]=2;
355                     } //reflection question
356                      
357                     //        Absorption loss 
358                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
359                       //printf("Got in absorption\n");
360                       ckovData[13]=20;
361                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
362                         ckovData[13]=11;
363                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
364                         ckovData[13]=12;
365                       if (gMC->CurrentVolID(copy) == gMC->VolId("META")) 
366                         ckovData[13]=13;
367                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP ")) 
368                         ckovData[13]=13;
369                       
370                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC")) 
371                         ckovData[13]=15;
372                       
373                       //        CsI inefficiency 
374                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
375                         ckovData[13]=16;
376                       }
377                       gMC->StopTrack();
378                       AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
379                     } //absorption question 
380                     
381                     
382                     //        Photon goes out of tracking scope 
383                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
384                       ckovData[13]=21;
385                       gMC->StopTrack();
386                       AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
387                     }   // energy treshold question         
388                   }  //number of mechanisms cycle
389                   /**********************End of evaluation************************/
390                 } //freon production question
391             } //energy interval question
392     } //cerenkov photon question
393       
394     /**************************************End of Production Parameters Storing*********************/ 
395     
396     
397     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/ 
398     
399     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
400       
401       if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
402         {
403           if (gMC->Edep() > 0.){
404                 gMC->TrackPosition(position);
405                 gMC->TrackMomentum(momentum);
406                 pos[0]=position(0);
407                 pos[1]=position(1);
408                 pos[2]=position(2);
409                 mom[0]=momentum(0);
410                 mom[1]=momentum(1);
411                 mom[2]=momentum(2);
412                 mom[3]=momentum(3);
413                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
414                 Double_t rt = TMath::Sqrt(tc);
415                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
416                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
417                 
418                 gMC->CurrentVolOffID(2,copy);
419                 vol[0]=copy;
420                 idvol=vol[0]-1;
421                 
422
423                 gMC->Gmtod(pos,localPos,1);
424
425                                                                     
426                 gMC->Gmtod(mom,localMom,2);
427
428                 
429                 gMC->CurrentVolOffID(2,copy);
430                 vol[0]=copy;
431                 idvol=vol[0]-1;
432
433                 
434                 ((AliRICHChamber*)fChambers->At(idvol))
435                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
436                 if(idvol<kNCH) {        
437                     ckovData[0] = gMC->TrackPid();        // particle type
438                     ckovData[1] = pos[0];                 // X-position for hit
439                     ckovData[2] = pos[1];                 // Y-position for hit
440                     ckovData[3] = pos[2];                 // Z-position for hit
441                     ckovData[4] = theta;                      // theta angle of incidence
442                     ckovData[5] = phi;                      // phi angle of incidence 
443                     ckovData[8] = (Float_t) fNsdigits;      // first sdigit
444                     ckovData[9] = -1;                       // last pad hit
445                     ckovData[13] = 4;                       // photon was detected
446                     ckovData[14] = mom[0];
447                     ckovData[15] = mom[1];
448                     ckovData[16] = mom[2];
449                     
450                     destep = gMC->Edep();
451                     gMC->SetMaxStep(kBig);
452                     cherenkovLoss  += destep;
453                     ckovData[7]=cherenkovLoss;
454                     
455                     //nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);//for photons in CsI kir
456                                     
457                     if (fNsdigits > (Int_t)ckovData[8]) {
458                         ckovData[8]= ckovData[8]+1;
459                         ckovData[9]= (Float_t) fNsdigits;
460                     }
461
462
463                     ckovData[17] = nPads;
464                     AliRICHhit *mipHit =  (AliRICHhit*) (fHits->UncheckedAt(0));
465                     if (mipHit)
466                       {
467                         mom[0] = current->Px();
468                         mom[1] = current->Py();
469                         mom[2] = current->Pz();
470                         Float_t mipPx = mipHit->MomX();
471                         Float_t mipPy = mipHit->MomY();
472                         Float_t mipPz = mipHit->MomZ();
473                         
474                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
475                         Float_t rt = TMath::Sqrt(r);
476                         Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz; 
477                         Float_t mipRt = TMath::Sqrt(mipR);
478                         if ((rt*mipRt) > 0)
479                           {
480                             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
481                           }
482                         else
483                           {
484                             coscerenkov = 0;
485                           }
486                         Float_t cherenkov = TMath::ACos(coscerenkov);
487                         ckovData[18]=cherenkov;
488                       }
489                     //if (sector != -1)
490                     //{
491                     AddHit(gAlice->GetCurrentTrackNumber(),vol,ckovData);
492                     AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
493                     //printf("Added One (5)!\n");
494                     //}
495                 }
496             }
497         }
498     }
499     
500     /***********************************************End of photon hits*********************************************/
501     
502
503     /**********************************************Charged particles treatment*************************************/
504
505     else if (gMC->TrackCharge()){//is MIP?
506         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
507           {
508             gMC->TrackMomentum(momentum);
509             mom[0]=momentum(0);
510             mom[1]=momentum(1);
511             mom[2]=momentum(2);
512             mom[3]=momentum(3);
513             hits [19] = mom[0];
514             hits [20] = mom[1];
515             hits [21] = mom[2];
516             fFreonProd=1;
517           }
518
519         if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)) {//is in GAP?
520 // Get current particle id (ipart), track position (pos)  and momentum (mom)
521             
522             gMC->CurrentVolOffID(3,copy);
523             vol[0]=copy;
524             idvol=vol[0]-1;
525             gMC->TrackPosition(position);
526             gMC->TrackMomentum(momentum);
527             pos[0]=position(0);
528             pos[1]=position(1);
529             pos[2]=position(2);
530             mom[0]=momentum(0);
531             mom[1]=momentum(1);
532             mom[2]=momentum(2);
533             mom[3]=momentum(3);
534
535             gMC->Gmtod(pos,localPos,1);
536             gMC->Gmtod(mom,localMom,2);
537             ipart  = gMC->TrackPid();
538             destep = gMC->Edep();step   = gMC->TrackStep();// momentum loss and steplength in last step
539             if(gMC->IsTrackEntering()){     // record hits when track enters ...
540
541                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
542                 Double_t rt = TMath::Sqrt(tc);
543                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
544                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
545                 
546
547                 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
548                 Double_t localRt = TMath::Sqrt(localTc);
549                 localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;                       
550                 localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;    
551                 
552                 hits[0] = Float_t(ipart);         // particle type
553                 hits[1] = localPos[0];                 // X-position for hit
554                 hits[2] = localPos[1];                 // Y-position for hit
555                 hits[3] = localPos[2];                 // Z-position for hit
556                 hits[4] = localTheta;                  // theta angle of incidence
557                 hits[5] = localPhi;                    // phi angle of incidence 
558                 hits[8] = (Float_t) fNsdigits;    // first sdigit
559                 hits[9] = -1;                     // last pad hit
560                 hits[13] = fFreonProd;           // did id hit the freon?
561                 hits[14] = mom[0];
562                 hits[15] = mom[1];
563                 hits[16] = mom[2];
564                 hits[18] = 0;               // dummy cerenkov angle
565
566                 tlength = 0;
567                 eloss   = 0;
568                 fFreonProd = 0;
569         
570                 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
571            
572                 
573                 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
574                 xhit    = localPos[0];
575                 yhit    = localPos[2];
576                 // Only if not trigger chamber
577                 if(idvol<kNCH) {
578                     //  Initialize hit position (cursor) in the segmentation model 
579                     ((AliRICHChamber*)fChambers->At(idvol))
580                         ->SigGenInit(localPos[0], localPos[2], localPos[1]);
581                 }
582             }
583             
584             // Calculate the charge induced on a pad (disintegration) in case 
585             // Mip left chamber ...
586             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
587                 gMC->SetMaxStep(kBig);
588                 eloss   += destep;
589                 tlength += step;
590                 
591                                 
592                 // Only if not trigger chamber
593                 if(idvol<kNCH) {
594                   if (eloss > 0) 
595                     {
596                       if(gMC->TrackPid() == kNeutron)
597                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
598                       //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); //for MIP kir
599                       hits[17] = nPads;
600                     }
601                 }
602                 
603                 hits[6]=tlength;
604                 hits[7]=eloss;
605                 if (fNsdigits > (Int_t)hits[8]) {
606                     hits[8]= hits[8]+1;
607                     hits[9]= (Float_t) fNsdigits;
608                 }
609                 
610                 new(lhits[fNhits++]) AliRICHhit(fIshunt,gAlice->GetCurrentTrackNumber(),vol,hits);
611                 eloss = 0; 
612                 // Check additional signal generation conditions 
613                 // defined by the segmentation
614                 // model (boundary crossing conditions) 
615             } else if 
616                 (((AliRICHChamber*)fChambers->At(idvol))
617                  ->SigGenCond(localPos[0], localPos[2], localPos[1]))
618             {
619                 ((AliRICHChamber*)fChambers->At(idvol))
620                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
621                 if (eloss > 0) 
622                   {
623                     if(gMC->TrackPid() == kNeutron)
624                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
625                     //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);//for N kir
626                     hits[17] = nPads;
627                     //printf("Npads:%d",NPads);
628                   }
629                 xhit     = localPos[0];
630                 yhit     = localPos[2]; 
631                 eloss    = destep;
632                 tlength += step ;
633                 //
634                 // nothing special  happened, add up energy loss
635             } else {        
636                 eloss   += destep;
637                 tlength += step ;
638             }
639         }//is in GAP?
640       }//is MIP?
641 }//void AliRICHv1::StepManager()