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 /*
  17   $Log$
  18   Revision 1.9  2000/05/31 08:19:38  jbarbosa
  19   Fixed bug in StepManager
  20
  21   Revision 1.8  2000/05/26 17:30:08  jbarbosa
  22   Cerenkov angle now stored within cerenkov data structure.
  23
  24   Revision 1.7  2000/05/18 10:31:36  jbarbosa
  25   Fixed positioning of spacers inside freon.
  26   Fixed positioning of proximity gap
  27   inside methane.
  28   Fixed cut on neutral particles in the StepManager.
  29
  30   Revision 1.6  2000/04/28 11:51:58  morsch
  31    Dimensions of arrays hits and Ckov_data corrected.
  32
  33   Revision 1.5  2000/04/19 13:28:46  morsch
  34   Major changes in geometry (parametrised), materials (updated) and
  35   step manager (diagnostics) (JB, AM)
  36
  37 */
  38
  39
  40
  41 //////////////////////////////////////////////////////////
  42 //  Manager and hits classes for set: RICH full version //
  43 //////////////////////////////////////////////////////////
  44
  45 #include <TTUBE.h>
  46 #include <TNode.h>
  47 #include <TRandom.h>
  48
  49 #include "AliRICHv1.h"
  50 #include "AliRun.h"
  51 #include "AliMC.h"
  52 #include "iostream.h"
  53 #include "AliCallf77.h"
  54 #include "AliConst.h"
  55 #include "AliPDG.h"
  56 #include "TGeant3.h"
  57
  58 ClassImp(AliRICHv1)
  59
  60 //___________________________________________
  61 AliRICHv1::AliRICHv1() : AliRICHv0()
  62 {
  63     //fChambers = 0;
  64 }
  65
  66 //___________________________________________
  67 AliRICHv1::AliRICHv1(const char *name, const char *title)
  68     : AliRICHv0(name,title)
  69 {
  70     fCkov_number=0;
  71     fFreon_prod=0;
  72
  73     fChambers = new TObjArray(7);
  74     for (Int_t i=0; i<7; i++) {
  75
  76         (*fChambers)[i] = new AliRICHChamber();
  77
  78     }
  79 }
  80
  81
  82 //___________________________________________
  83 void AliRICHv1::StepManager()
  84 {
  85     Int_t          copy, id;
  86     static Int_t   idvol;
  87     static Int_t   vol[2];
  88     Int_t          ipart;
  89     static Float_t hits[18];
  90     static Float_t Ckov_data[19];
  91     TLorentzVector Position;
  92     TLorentzVector Momentum;
  93     Float_t        pos[3];
  94     Float_t        mom[4];
  95     Float_t        Localpos[3];
  96     Float_t        Localmom[4];
  97     Float_t        Localtheta,Localphi;
  98     Float_t        theta,phi;
  99     Float_t        destep, step;
 100     Float_t        ranf[2];
 101     Int_t          NPads;
 102     Float_t        coscerenkov;
 103     static Float_t eloss, xhit, yhit, tlength;
 104     const  Float_t big=1.e10;
 105
 106     TClonesArray &lhits = *fHits;
 107     TGeant3 *geant3 = (TGeant3*) gMC;
 108     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
 109
 110  //if (current->Energy()>1)
 111    //{
 112
 113     // Only gas gap inside chamber
 114     // Tag chambers and record hits when track enters
 115
 116     idvol=-1;
 117     id=gMC->CurrentVolID(copy);
 118     Float_t cherenkov_loss=0;
 119     //gAlice->KeepTrack(gAlice->CurrentTrack());
 120
 121     gMC->TrackPosition(Position);
 122     pos[0]=Position(0);
 123     pos[1]=Position(1);
 124     pos[2]=Position(2);
 125     Ckov_data[1] = pos[0];                 // X-position for hit
 126     Ckov_data[2] = pos[1];                 // Y-position for hit
 127     Ckov_data[3] = pos[2];                 // Z-position for hit
 128     //Ckov_data[11] = gAlice->CurrentTrack();
 129
 130     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
 131
 132     /********************Store production parameters for Cerenkov photons************************/
 133 //is it a Cerenkov photon?
 134     if (gMC->TrackPid() == 50000050) {
 135
 136       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
 137         //{
 138           Float_t Ckov_energy = current->Energy();
 139           //energy interval for tracking
 140           if  (Ckov_energy > 5.6e-09 && Ckov_energy < 7.8e-09 )
 141             //if (Ckov_energy > 0)
 142             {
 143               if (gMC->IsTrackEntering()){                                     //is track entering?
 144                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 145                   {                                                          //is it in freo?
 146                     if (geant3->Gctrak()->nstep<1){                          //is it the first step?
 147                       Int_t mother = current->GetFirstMother();
 148
 149                       //printf("Second Mother:%d\n",current->GetSecondMother());
 150
 151                       Ckov_data[10] = mother;
 152                       Ckov_data[11] = gAlice->CurrentTrack();
 153                       Ckov_data[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
 154                       fCkov_number++;
 155                       fFreon_prod=1;
 156                       //printf("Index: %d\n",fCkov_number);
 157                     }    //first step question
 158                   }        //freo question
 159
 160                 if (geant3->Gctrak()->nstep<1){                                  //is it first step?
 161                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
 162                     {
 163                       Ckov_data[12] = 2;
 164                     }    //quarz question
 165                 }        //first step question
 166
 167                 //printf("Before %d\n",fFreon_prod);
 168               }   //track entering question
 169
 170               if (Ckov_data[12] == 1)                                        //was it produced in Freon?
 171                 //if (fFreon_prod == 1)
 172                 {
 173                   if (gMC->IsTrackEntering()){                                     //is track entering?
 174                     //printf("Got in");
 175                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
 176                       {
 177                         //printf("Got in\n");
 178                         gMC->TrackMomentum(Momentum);
 179                         mom[0]=Momentum(0);
 180                         mom[1]=Momentum(1);
 181                         mom[2]=Momentum(2);
 182                         mom[3]=Momentum(3);
 183                         // Z-position for hit
 184
 185
 186                         /**************** Photons lost in second grid have to be calculated by hand************/
 187
 188                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
 189                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
 190                         gMC->Rndm(ranf, 1);
 191                         //printf("grid calculation:%f\n",t);
 192                         if (ranf[0] > t) {
 193                           //geant3->StopTrack();
 194                           Ckov_data[13] = 5;
 195                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 196                           //printf("Lost one in grid\n");
 197                         }
 198                         /**********************************************************************************/
 199                       }    //gap
 200
 201                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
 202                       {
 203                         gMC->TrackMomentum(Momentum);
 204                         mom[0]=Momentum(0);
 205                         mom[1]=Momentum(1);
 206                         mom[2]=Momentum(2);
 207                         mom[3]=Momentum(3);
 208
 209                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
 210                         /***********************Cerenkov phtons (always polarised)*************************/
 211
 212                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
 213                         Float_t t = Fresnel(Ckov_energy*1e9,cophi,1);
 214                         gMC->Rndm(ranf, 1);
 215                         if (ranf[0] < t) {
 216                           //geant3->StopTrack();
 217                           Ckov_data[13] = 6;
 218                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 219                           //printf("Lost by Fresnel\n");
 220                         }
 221                         /**********************************************************************************/
 222                       }
 223                   } //track entering?
 224
 225
 226                   /********************Evaluation of losses************************/
 227                   /******************still in the old fashion**********************/
 228
 229                   Int_t i1 = geant3->Gctrak()->nmec;            //number of physics mechanisms acting on the particle
 230                   for (Int_t i = 0; i < i1; ++i) {
 231                     //        Reflection loss
 232                     if (geant3->Gctrak()->lmec[i] == 106) {        //was it reflected
 233                       Ckov_data[13]=10;
 234                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 235                         Ckov_data[13]=1;
 236                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 237                         Ckov_data[13]=2;
 238                       //geant3->StopTrack();
 239                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 240                     } //reflection question
 241
 242
 243                     //        Absorption loss
 244                     else if (geant3->Gctrak()->lmec[i] == 101) {              //was it absorbed?
 245                       Ckov_data[13]=20;
 246                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 247                         Ckov_data[13]=11;
 248                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 249                         Ckov_data[13]=12;
 250                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
 251                         Ckov_data[13]=13;
 252                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
 253                         Ckov_data[13]=13;
 254
 255                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
 256                         Ckov_data[13]=15;
 257
 258                       //        CsI inefficiency
 259                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
 260                         Ckov_data[13]=16;
 261                       }
 262                       //geant3->StopTrack();
 263                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 264                       //printf("Added cerenkov %d\n",fCkov_number);
 265                     } //absorption question
 266
 267
 268                     //        Photon goes out of tracking scope
 269                     else if (geant3->Gctrak()->lmec[i] == 30) {                 //is it below energy treshold?
 270                       Ckov_data[13]=21;
 271                       //geant3->StopTrack();
 272                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 273                     }   // energy treshold question
 274                   }  //number of mechanisms cycle
 275                   /**********************End of evaluation************************/
 276                 } //freon production question
 277             } //energy interval question
 278         //}//inside the proximity gap question
 279     } //cerenkov photon question
 280
 281     /**************************************End of Production Parameters Storing*********************/
 282
 283
 284     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
 285
 286     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
 287       //printf("Cerenkov\n");
 288         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
 289         {
 290
 291           if (gMC->Edep() > 0.){
 292                 gMC->TrackPosition(Position);
 293                 gMC->TrackMomentum(Momentum);
 294                 pos[0]=Position(0);
 295                 pos[1]=Position(1);
 296                 pos[2]=Position(2);
 297                 mom[0]=Momentum(0);
 298                 mom[1]=Momentum(1);
 299                 mom[2]=Momentum(2);
 300                 mom[3]=Momentum(3);
 301                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 302                 Double_t rt = TMath::Sqrt(tc);
 303                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 304                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 305                 gMC->Gmtod(pos,Localpos,1);
 306                 gMC->Gmtod(mom,Localmom,2);
 307
 308                 gMC->CurrentVolOffID(2,copy);
 309                 vol[0]=copy;
 310                 idvol=vol[0]-1;
 311
 312                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 313                         //->Sector(Localpos[0], Localpos[2]);
 314                 //printf("Sector:%d\n",sector);
 315
 316                 /*if (gMC->TrackPid() == 50000051){
 317                   fFeedbacks++;
 318                   printf("Feedbacks:%d\n",fFeedbacks);
 319                 }*/
 320
 321                 ((AliRICHChamber*) (*fChambers)[idvol])
 322                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 323                 if(idvol<7) {
 324                     Ckov_data[0] = gMC->TrackPid();        // particle type
 325                     Ckov_data[1] = pos[0];                 // X-position for hit
 326                     Ckov_data[2] = pos[1];                 // Y-position for hit
 327                     Ckov_data[3] = pos[2];                 // Z-position for hit
 328                     Ckov_data[4] = theta;                      // theta angle of incidence
 329                     Ckov_data[5] = phi;                      // phi angle of incidence
 330                     Ckov_data[8] = (Float_t) fNPadHits;      // first padhit
 331                     Ckov_data[9] = -1;                       // last pad hit
 332                     Ckov_data[13] = 4;                       // photon was detected
 333                     Ckov_data[14] = mom[0];
 334                     Ckov_data[15] = mom[1];
 335                     Ckov_data[16] = mom[2];
 336
 337                     destep = gMC->Edep();
 338                     gMC->SetMaxStep(big);
 339                     cherenkov_loss  += destep;
 340                     Ckov_data[7]=cherenkov_loss;
 341
 342                     NPads = MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
 343                     if (fNPadHits > (Int_t)Ckov_data[8]) {
 344                         Ckov_data[8]= Ckov_data[8]+1;
 345                         Ckov_data[9]= (Float_t) fNPadHits;
 346                     }
 347
 348                     Ckov_data[17] = NPads;
 349                     //printf("Npads:%d",NPads);
 350
 351                     //TClonesArray *Hits = RICH->Hits();
 352                     AliRICHHit *mipHit =  (AliRICHHit*) (fHits->UncheckedAt(0));
 353                     if (mipHit)
 354                       {
 355                         mom[0] = current->Px();
 356                         mom[1] = current->Py();
 357                         mom[2] = current->Pz();
 358                         Float_t Mip_px = mipHit->fMomX;
 359                         Float_t Mip_py = mipHit->fMomY;
 360                         Float_t Mip_pz = mipHit->fMomZ;
 361
 362                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
 363                         Float_t rt = TMath::Sqrt(r);
 364                         Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
 365                         Float_t Mip_rt = TMath::Sqrt(Mip_r);
 366                         if ((rt*Mip_rt) > 0)
 367                           {
 368                             coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt);
 369                           }
 370                         else
 371                           {
 372                             coscerenkov = 0;
 373                           }
 374                         Float_t cherenkov = TMath::ACos(coscerenkov);
 375                         Ckov_data[18]=cherenkov;
 376                       }
 377                     //if (sector != -1)
 378                     //{
 379                     AddHit(gAlice->CurrentTrack(),vol,Ckov_data);
 380                     AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
 381                     //}
 382                 }
 383             }
 384         }
 385     }
 386
 387     /***********************************************End of photon hits*********************************************/
 388
 389
 390     /**********************************************Charged particles treatment*************************************/
 391
 392     else if (gMC->TrackCharge())
 393     //else if (1 == 1)
 394       {
 395 //If MIP
 396         /*if (gMC->IsTrackEntering())
 397           {
 398             hits[13]=20;//is track entering?
 399           }*/
 400         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 401           {
 402             fFreon_prod=1;
 403           }
 404
 405         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
 406 // Get current particle id (ipart), track position (pos)  and momentum (mom)
 407
 408             gMC->CurrentVolOffID(3,copy);
 409             vol[0]=copy;
 410             idvol=vol[0]-1;
 411
 412             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 413                         //->Sector(Localpos[0], Localpos[2]);
 414             //printf("Sector:%d\n",sector);
 415
 416             gMC->TrackPosition(Position);
 417             gMC->TrackMomentum(Momentum);
 418             pos[0]=Position(0);
 419             pos[1]=Position(1);
 420             pos[2]=Position(2);
 421             mom[0]=Momentum(0);
 422             mom[1]=Momentum(1);
 423             mom[2]=Momentum(2);
 424             mom[3]=Momentum(3);
 425             gMC->Gmtod(pos,Localpos,1);
 426             gMC->Gmtod(mom,Localmom,2);
 427
 428             ipart  = gMC->TrackPid();
 429             //
 430             // momentum loss and steplength in last step
 431             destep = gMC->Edep();
 432             step   = gMC->TrackStep();
 433
 434             //
 435             // record hits when track enters ...
 436             if( gMC->IsTrackEntering()) {
 437 //              gMC->SetMaxStep(fMaxStepGas);
 438                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 439                 Double_t rt = TMath::Sqrt(tc);
 440                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 441                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 442
 443
 444                 Double_t Localtc = Localmom[0]*Localmom[0]+Localmom[2]*Localmom[2];
 445                 Double_t Localrt = TMath::Sqrt(Localtc);
 446                 Localtheta   = Float_t(TMath::ATan2(Localrt,Double_t(Localmom[1])))*kRaddeg;
 447                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[2]),Double_t(Localmom[0])))*kRaddeg;
 448
 449                 hits[0] = Float_t(ipart);         // particle type
 450                 hits[1] = Localpos[0];                 // X-position for hit
 451                 hits[2] = Localpos[1];                 // Y-position for hit
 452                 hits[3] = Localpos[2];                 // Z-position for hit
 453                 hits[4] = Localtheta;                  // theta angle of incidence
 454                 hits[5] = Localphi;                    // phi angle of incidence
 455                 hits[8] = (Float_t) fNPadHits;    // first padhit
 456                 hits[9] = -1;                     // last pad hit
 457                 hits[13] = fFreon_prod;           // did id hit the freon?
 458                 hits[14] = mom[0];
 459                 hits[15] = mom[1];
 460                 hits[16] = mom[2];
 461
 462                 tlength = 0;
 463                 eloss   = 0;
 464                 fFreon_prod = 0;
 465
 466                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
 467
 468
 469                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
 470                 xhit    = Localpos[0];
 471                 yhit    = Localpos[2];
 472                 // Only if not trigger chamber
 473                 if(idvol<7) {
 474                     //
 475                     //  Initialize hit position (cursor) in the segmentation model
 476                     ((AliRICHChamber*) (*fChambers)[idvol])
 477                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 478                 }
 479             }
 480
 481             //
 482             // Calculate the charge induced on a pad (disintegration) in case
 483             //
 484             // Mip left chamber ...
 485             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 486                 gMC->SetMaxStep(big);
 487                 eloss   += destep;
 488                 tlength += step;
 489
 490
 491                 // Only if not trigger chamber
 492                 if(idvol<7) {
 493                   if (eloss > 0)
 494                     {
 495                       if(gMC->TrackPid() == kNeutron)
 496                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 497                       NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
 498                       hits[17] = NPads;
 499                       //printf("Npads:%d",NPads);
 500                     }
 501                 }
 502
 503                 hits[6]=tlength;
 504                 hits[7]=eloss;
 505                 if (fNPadHits > (Int_t)hits[8]) {
 506                     hits[8]= hits[8]+1;
 507                     hits[9]= (Float_t) fNPadHits;
 508                 }
 509
 510                 //if(sector !=-1)
 511                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 512                 eloss = 0;
 513                 //
 514                 // Check additional signal generation conditions
 515                 // defined by the segmentation
 516                 // model (boundary crossing conditions)
 517             } else if
 518                 (((AliRICHChamber*) (*fChambers)[idvol])
 519                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
 520             {
 521                 ((AliRICHChamber*) (*fChambers)[idvol])
 522                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 523                 if (eloss > 0)
 524                   {
 525                     if(gMC->TrackPid() == kNeutron)
 526                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 527                     NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
 528                     hits[17] = NPads;
 529                     //printf("Npads:%d",NPads);
 530                   }
 531                 xhit     = Localpos[0];
 532                 yhit     = Localpos[2];
 533                 eloss    = destep;
 534                 tlength += step ;
 535                 //
 536                 // nothing special  happened, add up energy loss
 537             } else {
 538                 eloss   += destep;
 539                 tlength += step ;
 540             }
 541         }
 542       }
 543     /*************************************************End of MIP treatment**************************************/
 544    //}
 545 }