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