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