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