RICH/AliRICHv3.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 #include <Riostream.h>
  18
  19 #include <TBRIK.h>
  20 #include <TGeometry.h>
  21 #include <TLorentzVector.h>
  22 #include <TNode.h>
  23 #include <TParticle.h>
  24 #include <TVector3.h>
  25 #include <TVirtualMC.h>
  26 #include <TPDGCode.h> //for kNuetron
  27 #include <TCanvas.h>
  28 #include <TF1.h>
  29 #include <TH1.h>
  30 #include <TH2.h>
  31 #include <TStyle.h>
  32
  33 #include "AliConst.h"
  34 #include "AliMagF.h"
  35 #include "AliPDG.h"
  36 #include "AliRICHGeometry.h"
  37 #include "AliRICHResponse.h"
  38 #include "AliRICHSegmentationV1.h"
  39 #include "AliRICHv3.h"
  40 #include "AliRun.h"
  41 #include "AliMC.h"
  42
  43
  44 ClassImp(AliRICHv3)
  45
  46 //______________________________________________________________
  47 //    Implementation of the RICH version 3 with azimuthal rotation
  48
  49
  50 AliRICHv3::AliRICHv3(const char *sName, const char *sTitle)
  51           :AliRICH(sName,sTitle)
  52 {
  53 // The named ctor currently creates a single copy of
  54 // AliRICHGeometry AliRICHSegmentationV1 AliRICHResponse
  55 // and initialises the corresponding models of all 7 chambers with these stuctures.
  56 // Note: all chambers share the single copy of models. MUST be changed later (???).
  57   if(GetDebug())Info("named ctor","Start.");
  58
  59    fCkovNumber=fFreonProd=0;
  60
  61    AliRICHGeometry       *pRICHGeometry    =new AliRICHGeometry;           // ??? to be moved to AlRICHChamber::named ctor
  62    AliRICHSegmentationV1 *pRICHSegmentation=new AliRICHSegmentationV1;     // ??? to be moved to AlRICHChamber::named ctor
  63    AliRICHResponse       *pRICHResponse    =new AliRICHResponse;           // ??? to be moved to AlRICHChamber::named ctor
  64
  65    for (Int_t i=1; i<=kNCH; i++){
  66       SetGeometryModel(i,pRICHGeometry);
  67       SetSegmentationModel(i,pRICHSegmentation);
  68       SetResponseModel(i,pRICHResponse);
  69       C(i)->Init(i); // ??? to be removed
  70    }
  71   if(GetDebug())Info("named ctor","Stop.");
  72 }//AliRICHv3::ctor(const char *pcName, const char *pcTitle)
  73
  74 AliRICHv3::~AliRICHv3()
  75 {
  76 // Dtor deletes RICH models. In future (???) AliRICHChamber will be responsible for that.
  77    if(GetDebug()) cout<<ClassName()<<"::dtor()>\n";
  78
  79    if(fChambers) {
  80      AliRICHChamber *ch =C(1);
  81      if(ch) {
  82        delete ch->GetGeometryModel();
  83        delete ch->GetResponseModel();
  84        delete ch->GetSegmentationModel();
  85      }
  86    }
  87 }//AliRICHv3::dtor()
  88 //______________________________________________________________________________
  89 void AliRICHv3::StepManager()
  90 {//Full Step Manager
  91
  92     Int_t          copy, id;
  93     static Int_t   idvol;
  94     static Int_t   vol[2];
  95     Int_t          ipart;
  96     static Float_t hits[22];
  97     static Float_t ckovData[19];
  98     TLorentzVector position;
  99     TLorentzVector momentum;
 100     Float_t        pos[3];
 101     Float_t        mom[4];
 102     Float_t        localPos[3];
 103     Float_t        localMom[4];
 104     Float_t        localTheta,localPhi;
 105     Float_t        theta,phi;
 106     Float_t        destep, step;
 107     Double_t        ranf[2];
 108     Float_t        coscerenkov;
 109     static Float_t eloss, xhit, yhit, tlength;
 110     const  Float_t kBig=1.e10;
 111
 112     TClonesArray &lhits = *fHits;
 113     TParticle *current = (TParticle*)(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
 114
 115  //if (current->Energy()>1)
 116    //{
 117
 118     // Only gas gap inside chamber
 119     // Tag chambers and record hits when track enters
 120
 121
 122     id=gMC->CurrentVolID(copy);
 123     idvol = copy-1;
 124     Float_t cherenkovLoss=0;
 125     //gAlice->KeepTrack(gAlice->GetCurrentTrackNumber());
 126
 127     gMC->TrackPosition(position);
 128     pos[0]=position(0);
 129     pos[1]=position(1);
 130     pos[2]=position(2);
 131     //bzero((char *)ckovData,sizeof(ckovData)*19);
 132     ckovData[1] = pos[0];                 // X-position for hit
 133     ckovData[2] = pos[1];                 // Y-position for hit
 134     ckovData[3] = pos[2];                 // Z-position for hit
 135     ckovData[6] = 0;                      // dummy track length
 136     //ckovData[11] = gAlice->GetCurrentTrackNumber();
 137
 138     //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->GetCurrentTrackNumber());
 139
 140     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
 141
 142     /********************Store production parameters for Cerenkov photons************************/
 143 //is it a Cerenkov photon?
 144     if (gMC->TrackPid() == 50000050) {
 145
 146       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
 147         //{
 148           Float_t ckovEnergy = current->Energy();
 149           //energy interval for tracking
 150           if  (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
 151             //if (ckovEnergy > 0)
 152             {
 153               if (gMC->IsTrackEntering()){        //is track entering?
 154                 //printf("Track entered (1)\n");
 155                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 156                   {                                                          //is it in freo?
 157                     if (gMC->IsNewTrack()){                          //is it the first step?
 158                       //printf("I'm in!\n");
 159                       Int_t mother = current->GetFirstMother();
 160
 161                       //printf("Second Mother:%d\n",current->GetSecondMother());
 162
 163                       ckovData[10] = mother;
 164                       ckovData[11] = gAlice->GetMCApp()->GetCurrentTrackNumber();
 165                       ckovData[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
 166                       //printf("Produced in FREO\n");
 167                       fCkovNumber++;
 168                       fFreonProd=1;
 169                       //printf("Index: %d\n",fCkovNumber);
 170                     }    //first step question
 171                   }        //freo question
 172
 173                 if (gMC->IsNewTrack()){                                  //is it first step?
 174                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
 175                     {
 176                       ckovData[12] = 2;
 177                       //printf("Produced in QUAR\n");
 178                     }    //quarz question
 179                 }        //first step question
 180
 181                 //printf("Before %d\n",fFreonProd);
 182               }   //track entering question
 183
 184               if (ckovData[12] == 1)                                        //was it produced in Freon?
 185                 //if (fFreonProd == 1)
 186                 {
 187                   if (gMC->IsTrackEntering()){                                     //is track entering?
 188                     //printf("Track entered (2)\n");
 189                     //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
 190                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
 191                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
 192                       {
 193                         //printf("Got in META\n");
 194                         gMC->TrackMomentum(momentum);
 195                         mom[0]=momentum(0);
 196                         mom[1]=momentum(1);
 197                         mom[2]=momentum(2);
 198                         mom[3]=momentum(3);
 199
 200                         gMC->Gmtod(mom,localMom,2);
 201                         Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
 202                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
 203                         /**************** Photons lost in second grid have to be calculated by hand************/
 204                         gMC->GetRandom()->RndmArray(1,ranf);
 205                         if (ranf[0] > t) {
 206                           gMC->StopTrack();
 207                           ckovData[13] = 5;
 208                           AddCerenkov(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 209                           //printf("Added One (1)!\n");
 210                           //printf("Lost one in grid\n");
 211                         }
 212                         /**********************************************************************************/
 213                       }    //gap
 214
 215                     //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
 216                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
 217                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
 218                       {
 219                         //printf("Got in CSI\n");
 220                         gMC->TrackMomentum(momentum);
 221                         mom[0]=momentum(0);
 222                         mom[1]=momentum(1);
 223                         mom[2]=momentum(2);
 224                         mom[3]=momentum(3);
 225
 226                         gMC->Gmtod(mom,localMom,2);
 227                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
 228                         /***********************Cerenkov phtons (always polarised)*************************/
 229                         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 230                         Double_t localRt = TMath::Sqrt(localTc);
 231                         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
 232                         Double_t cotheta = TMath::Abs(cos(localTheta));
 233                         Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
 234                             gMC->GetRandom()->RndmArray(1,ranf);
 235                             if (ranf[0] < t) {
 236                               gMC->StopTrack();
 237                               ckovData[13] = 6;
 238                               AddCerenkov(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 239
 240                               //printf("Added One (2)!\n");
 241                               //printf("Lost by Fresnel\n");
 242                             }
 243                             /**********************************************************************************/
 244                       }
 245                   } //track entering?
 246
 247
 248                   /********************Evaluation of losses************************/
 249                   /******************still in the old fashion**********************/
 250
 251                   TArrayI procs;
 252                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
 253                   for (Int_t i = 0; i < i1; ++i) {
 254                     //        Reflection loss
 255                     if (procs[i] == kPLightReflection) {        //was it reflected
 256                       ckovData[13]=10;
 257                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 258                         ckovData[13]=1;
 259                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 260                         ckovData[13]=2;
 261                       //gMC->StopTrack();
 262                       //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 263                     } //reflection question
 264
 265                     //        Absorption loss
 266                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
 267                       //printf("Got in absorption\n");
 268                       ckovData[13]=20;
 269                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 270                         ckovData[13]=11;
 271                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 272                         ckovData[13]=12;
 273                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
 274                         ckovData[13]=13;
 275                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
 276                         ckovData[13]=13;
 277
 278                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
 279                         ckovData[13]=15;
 280
 281                       //        CsI inefficiency
 282                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
 283                         ckovData[13]=16;
 284                       }
 285                       gMC->StopTrack();
 286                       AddCerenkov(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 287                       //printf("Added One (3)!\n");
 288                       //printf("Added cerenkov %d\n",fCkovNumber);
 289                     } //absorption question
 290
 291
 292                     //        Photon goes out of tracking scope
 293                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
 294                       ckovData[13]=21;
 295                       gMC->StopTrack();
 296                       AddCerenkov(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 297                       //printf("Added One (4)!\n");
 298                     }   // energy treshold question
 299                   }  //number of mechanisms cycle
 300                   /**********************End of evaluation************************/
 301                 } //freon production question
 302             } //energy interval question
 303         //}//inside the proximity gap question
 304     } //cerenkov photon question
 305
 306     /**************************************End of Production Parameters Storing*********************/
 307
 308
 309     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
 310
 311     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
 312       //printf("Cerenkov\n");
 313
 314       //if (gMC->TrackPid() == 50000051)
 315         //printf("Tracking a feedback\n");
 316
 317       if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
 318         {
 319           //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
 320           //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
 321           //printf("Got in CSI\n");
 322           //printf("Tracking a %d\n",gMC->TrackPid());
 323           if (gMC->Edep() > 0.){
 324                 gMC->TrackPosition(position);
 325                 gMC->TrackMomentum(momentum);
 326                 pos[0]=position(0);
 327                 pos[1]=position(1);
 328                 pos[2]=position(2);
 329                 mom[0]=momentum(0);
 330                 mom[1]=momentum(1);
 331                 mom[2]=momentum(2);
 332                 mom[3]=momentum(3);
 333                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 334                 Double_t rt = TMath::Sqrt(tc);
 335                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 336                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 337
 338                 gMC->CurrentVolOffID(2,copy);
 339                 vol[0]=copy;
 340                 idvol=vol[0]-1;
 341
 342
 343                 gMC->Gmtod(pos,localPos,1);
 344
 345                 //Chamber(idvol).GlobaltoLocal(pos,localPos);
 346
 347                 gMC->Gmtod(mom,localMom,2);
 348
 349                 //Chamber(idvol).GlobaltoLocal(mom,localMom);
 350
 351                 gMC->CurrentVolOffID(2,copy);
 352                 vol[0]=copy;
 353                 idvol=vol[0]-1;
 354
 355                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 356                         //->Sector(localPos[0], localPos[2]);
 357                 //printf("Sector:%d\n",sector);
 358
 359                 /*if (gMC->TrackPid() == 50000051){
 360                   fFeedbacks++;
 361                   printf("Feedbacks:%d\n",fFeedbacks);
 362                 }*/
 363
 364         //PH            ((AliRICHChamber*) (*fChambers)[idvol])
 365                 ((AliRICHChamber*)fChambers->At(idvol))
 366                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
 367                 if(idvol<kNCH) {
 368                     ckovData[0] = gMC->TrackPid();        // particle type
 369                     ckovData[1] = pos[0];                 // X-position for hit
 370                     ckovData[2] = pos[1];                 // Y-position for hit
 371                     ckovData[3] = pos[2];                 // Z-position for hit
 372                     ckovData[4] = theta;                      // theta angle of incidence
 373                     ckovData[5] = phi;                      // phi angle of incidence
 374                     ckovData[8] = (Float_t) fNsdigits;      // first sdigit
 375                     ckovData[9] = -1;                       // last pad hit
 376                     ckovData[13] = 4;                       // photon was detected
 377                     ckovData[14] = mom[0];
 378                     ckovData[15] = mom[1];
 379                     ckovData[16] = mom[2];
 380
 381                     destep = gMC->Edep();
 382                     gMC->SetMaxStep(kBig);
 383                     cherenkovLoss  += destep;
 384                     ckovData[7]=cherenkovLoss;
 385
 386                     ckovData[17] = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kPhoton);//for photons in CsI
 387
 388                     if (fNsdigits > (Int_t)ckovData[8]) {
 389                         ckovData[8]= ckovData[8]+1;
 390                         ckovData[9]= (Float_t) fNsdigits;
 391                     }
 392
 393
 394                     //TClonesArray *Hits = RICH->Hits();
 395                     AliRICHhit *mipHit =  (AliRICHhit*) (fHits->UncheckedAt(0));
 396                     if (mipHit)
 397                       {
 398                         mom[0] = current->Px();
 399                         mom[1] = current->Py();
 400                         mom[2] = current->Pz();
 401                         Float_t mipPx = mipHit->MomX();
 402                         Float_t mipPy = mipHit->MomY();
 403                         Float_t mipPz = mipHit->MomZ();
 404
 405                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
 406                         Float_t rt = TMath::Sqrt(r);
 407                         Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
 408                         Float_t mipRt = TMath::Sqrt(mipR);
 409                         if ((rt*mipRt) > 0)
 410                           {
 411                             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
 412                           }
 413                         else
 414                           {
 415                             coscerenkov = 0;
 416                           }
 417                         Float_t cherenkov = TMath::ACos(coscerenkov);
 418                         ckovData[18]=cherenkov;
 419                       }
 420                     //if (sector != -1)
 421                     //{
 422                     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 423                     AddCerenkov(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);
 424                     //printf("Added One (5)!\n");
 425                     //}
 426                 }
 427             }
 428         }
 429     }
 430
 431     /***********************************************End of photon hits*********************************************/
 432
 433
 434     /**********************************************Charged particles treatment*************************************/
 435
 436     else if (gMC->TrackCharge()){
 437 //If MIP
 438         /*if (gMC->IsTrackEntering())
 439           {
 440             hits[13]=20;//is track entering?
 441           }*/
 442         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 443           {
 444             gMC->TrackMomentum(momentum);
 445             mom[0]=momentum(0);
 446             mom[1]=momentum(1);
 447             mom[2]=momentum(2);
 448             mom[3]=momentum(3);
 449             hits [19] = mom[0];
 450             hits [20] = mom[1];
 451             hits [21] = mom[2];
 452             fFreonProd=1;
 453           }
 454
 455         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {//is in GAP?
 456 // Get current particle id (ipart), track position (pos)  and momentum (mom)
 457
 458             gMC->CurrentVolOffID(3,copy);
 459             vol[0]=copy;
 460             idvol=vol[0]-1;
 461
 462             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 463                         //->Sector(localPos[0], localPos[2]);
 464             //printf("Sector:%d\n",sector);
 465
 466             gMC->TrackPosition(position);
 467             gMC->TrackMomentum(momentum);
 468             pos[0]=position(0);
 469             pos[1]=position(1);
 470             pos[2]=position(2);
 471             mom[0]=momentum(0);
 472             mom[1]=momentum(1);
 473             mom[2]=momentum(2);
 474             mom[3]=momentum(3);
 475
 476             gMC->Gmtod(pos,localPos,1);
 477
 478             //Chamber(idvol).GlobaltoLocal(pos,localPos);
 479
 480             gMC->Gmtod(mom,localMom,2);
 481
 482             //Chamber(idvol).GlobaltoLocal(mom,localMom);
 483
 484             ipart  = gMC->TrackPid();
 485             //
 486             // momentum loss and steplength in last step
 487             destep = gMC->Edep();
 488             step   = gMC->TrackStep();
 489
 490             //
 491             // record hits when track enters ...
 492             if( gMC->IsTrackEntering()) {
 493 //              gMC->SetMaxStep(fMaxStepGas);
 494                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 495                 Double_t rt = TMath::Sqrt(tc);
 496                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 497                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 498
 499
 500                 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 501                 Double_t localRt = TMath::Sqrt(localTc);
 502                 localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
 503                 localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
 504
 505                 hits[0] = Float_t(ipart);         // particle type
 506                 hits[1] = localPos[0];                 // X-position for hit
 507                 hits[2] = localPos[1];                 // Y-position for hit
 508                 hits[3] = localPos[2];                 // Z-position for hit
 509                 hits[4] = localTheta;                  // theta angle of incidence
 510                 hits[5] = localPhi;                    // phi angle of incidence
 511                 hits[8] = (Float_t) fNsdigits;    // first sdigit
 512                 hits[9] = -1;                     // last pad hit
 513                 hits[13] = fFreonProd;           // did id hit the freon?
 514                 hits[14] = mom[0];
 515                 hits[15] = mom[1];
 516                 hits[16] = mom[2];
 517                 hits[18] = 0;               // dummy cerenkov angle
 518
 519                 tlength = 0;
 520                 eloss   = 0;
 521                 fFreonProd = 0;
 522
 523                 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
 524
 525
 526                 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
 527                 xhit    = localPos[0];
 528                 yhit    = localPos[2];
 529                 // Only if not trigger chamber
 530                 if(idvol<kNCH) {
 531                     //
 532                     //  Initialize hit position (cursor) in the segmentation model
 533           //PH              ((AliRICHChamber*) (*fChambers)[idvol])
 534                     ((AliRICHChamber*)fChambers->At(idvol))
 535                         ->SigGenInit(localPos[0], localPos[2], localPos[1]);
 536                 }
 537             }
 538
 539             //
 540             // Calculate the charge induced on a pad (disintegration) in case
 541             //
 542             // Mip left chamber ...
 543             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 544                 gMC->SetMaxStep(kBig);
 545                 eloss   += destep;
 546                 tlength += step;
 547
 548
 549                 // Only if not trigger chamber
 550                 if(idvol<kNCH) {
 551                   if (eloss > 0)
 552                     {
 553                       if(gMC->TrackPid() == kNeutron)
 554                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 555                       hits[17] = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); //for MIP
 556                     }
 557                 }
 558
 559                 hits[6]=tlength;
 560                 hits[7]=eloss;
 561                 if (fNsdigits > (Int_t)hits[8]) {
 562                     hits[8]= hits[8]+1;
 563                     hits[9]= (Float_t) fNsdigits;
 564                 }
 565
 566                 //if(sector !=-1)
 567                 new(lhits[fNhits++]) AliRICHhit(fIshunt,gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);
 568                 eloss = 0;
 569                 //
 570                 // Check additional signal generation conditions
 571                 // defined by the segmentation
 572                 // model (boundary crossing conditions)
 573             }else if(((AliRICHChamber*)fChambers->At(idvol))->SigGenCond(localPos[0], localPos[2], localPos[1])){
 574                 ((AliRICHChamber*)fChambers->At(idvol))->SigGenInit(localPos[0], localPos[2], localPos[1]);
 575                 if (eloss > 0)
 576                   {
 577                     if(gMC->TrackPid() == kNeutron)
 578                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 579                     hits[17] = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);//for n
 580                   }
 581                 xhit     = localPos[0];
 582                 yhit     = localPos[2];
 583                 eloss    = destep;
 584                 tlength += step ;
 585                 //
 586                 // nothing special  happened, add up energy loss
 587             } else {
 588                 eloss   += destep;
 589                 tlength += step ;
 590             }
 591         }//is in GAP?
 592       }//is MIP?
 593     /*************************************************End of MIP treatment**************************************/
 594 }//void AliRICHv3::StepManager()
 595 //__________________________________________________________________________________________________
 596 Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
 597 {//calls the charge disintegration method of the current chamber and adds all generated sdigits to the list of digits
 598
 599    Float_t newclust[4][500];
 600    Int_t clhits[5];
 601    Int_t iNdigits;
 602    clhits[0]=fNhits+1;
 603
 604   ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, iNdigits,newclust, res);
 605
 606   for (Int_t i=0; i<iNdigits; i++) {
 607     if (Int_t(newclust[0][i]) > 0) {
 608             clhits[1] = Int_t(newclust[0][i]);//  Cluster Charge
 609             clhits[2] = Int_t(newclust[1][i]);//  Pad: ix
 610             clhits[3] = Int_t(newclust[2][i]);//  Pad: iy
 611             clhits[4] = Int_t(newclust[3][i]);//  Pad: chamber sector
 612             AddSpecialOld(clhits);
 613         }
 614     }
 615   return iNdigits;
 616 }//Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
 617 //__________________________________________________________________________________________________