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(0);
 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;
 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()