RICH/AliRICHv0.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.5  2000/04/19 13:28:46  morsch
  19   Major changes in geometry (parametrised), materials (updated) and
  20   step manager (diagnostics) (JB, AM)
  21
  22 */
  23
  24
  25
  26 ////////////////////////////////////////////////////////
  27 //  Manager and hits classes for set:RICH version 0    //
  28 /////////////////////////////////////////////////////////
  29
  30 #include <TTUBE.h>
  31 #include <TNode.h>
  32 #include <TRandom.h>
  33
  34 #include "AliRICHv0.h"
  35 #include "AliRun.h"
  36 #include "AliMC.h"
  37 #include "iostream.h"
  38 #include "AliCallf77.h"
  39 #include "AliConst.h"
  40 #include "AliPDG.h"
  41 #include "TGeant3.h"
  42
  43 ClassImp(AliRICHv0)
  44
  45 //___________________________________________
  46 AliRICHv0::AliRICHv0() : AliRICH()
  47 {
  48     //fChambers = 0;
  49 }
  50
  51 //___________________________________________
  52 AliRICHv0::AliRICHv0(const char *name, const char *title)
  53     : AliRICH(name,title)
  54 {
  55     fCkov_number=0;
  56     fFreon_prod=0;
  57
  58     fChambers = new TObjArray(7);
  59     for (Int_t i=0; i<7; i++) {
  60
  61         (*fChambers)[i] = new AliRICHChamber();
  62
  63     }
  64 }
  65
  66
  67 //___________________________________________
  68 void AliRICHv0::CreateGeometry()
  69 {
  70     //
  71     // Create the geometry for RICH version 1
  72     //
  73     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
  74     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
  75     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
  76     //
  77     //Begin_Html
  78     /*
  79       <img src="picts/AliRICHv1.gif">
  80     */
  81     //End_Html
  82     //Begin_Html
  83     /*
  84       <img src="picts/AliRICHv1Tree.gif">
  85     */
  86     //End_Html
  87
  88   AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
  89   AliRICHSegmentation*  segmentation;
  90   AliRICHGeometry*  geometry;
  91   AliRICHChamber*       iChamber;
  92
  93   iChamber = &(RICH->Chamber(0));
  94   segmentation=iChamber->GetSegmentationModel(0);
  95   geometry=iChamber->GetGeometryModel();
  96
  97
  98     Int_t *idtmed = fIdtmed->GetArray()-999;
  99
 100     Int_t i;
 101     Float_t zs;
 102     Int_t idrotm[1099];
 103     Float_t par[3];
 104
 105     // --- Define the RICH detector
 106     //     External aluminium box
 107     par[0] = 71.1;
 108     par[1] = 11.5;                 //Original Settings
 109     par[2] = 73.15;
 110     /*par[0] = 73.15;
 111     par[1] = 11.5;
 112     par[2] = 71.1;*/
 113     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 114
 115     //     Sensitive part of the whole RICH
 116     par[0] = 64.8;
 117     par[1] = 11.5;                 //Original Settings
 118     par[2] = 66.55;
 119     /*par[0] = 66.55;
 120     par[1] = 11.5;
 121     par[2] = 64.8;*/
 122     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
 123
 124     //     Honeycomb
 125     par[0] = 63.1;
 126     par[1] = .188;                 //Original Settings
 127     par[2] = 66.55;
 128     /*par[0] = 66.55;
 129     par[1] = .188;
 130     par[2] = 63.1;*/
 131     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
 132
 133     //     Aluminium sheet
 134     par[0] = 63.1;
 135     par[1] = .025;                 //Original Settings
 136     par[2] = 66.55;
 137     /*par[0] = 66.5;
 138     par[1] = .025;
 139     par[2] = 63.1;*/
 140     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
 141
 142     //     Quartz
 143     par[0] = geometry->GetQuartzWidth()/2;
 144     par[1] = geometry->GetQuartzThickness()/2;
 145     par[2] = geometry->GetQuartzLength()/2;
 146     /*par[0] = 63.1;
 147     par[1] = .25;                  //Original Settings
 148     par[2] = 65.5;*/
 149     /*par[0] = geometry->GetQuartzWidth()/2;
 150     par[1] = geometry->GetQuartzThickness()/2;
 151     par[2] = geometry->GetQuartzLength()/2;*/
 152     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
 153     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
 154
 155     //     Spacers (cylinders)
 156     par[0] = 0.;
 157     par[1] = .5;
 158     par[2] = geometry->GetFreonThickness()/2;
 159     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
 160
 161     //     Opaque quartz
 162     par[0] = 61.95;
 163     par[1] = .2;                   //Original Settings
 164     par[2] = 66.5;
 165     /*par[0] = 66.5;
 166     par[1] = .2;
 167     par[2] = 61.95;*/
 168     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
 169
 170     //     Frame of opaque quartz
 171     par[0] = geometry->GetOuterFreonWidth()/2;
 172     par[1] = geometry->GetFreonThickness()/2;
 173     par[2] = geometry->GetOuterFreonLength()/2 + 1;
 174     /*par[0] = 20.65;
 175     par[1] = .5;                   //Original Settings
 176     par[2] = 66.5;*/
 177     /*par[0] = 66.5;
 178     par[1] = .5;
 179     par[2] = 20.65;*/
 180     gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
 181
 182     par[0] = geometry->GetInnerFreonWidth()/2;
 183     par[1] = geometry->GetFreonThickness()/2;
 184     par[2] = geometry->GetInnerFreonLength()/2 + 1;
 185     gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
 186
 187     //     Little bar of opaque quartz
 188     par[0] = .275;
 189     par[1] = geometry->GetQuartzThickness()/2;
 190     par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
 191     /*par[0] = .275;
 192     par[1] = .25;                   //Original Settings
 193     par[2] = 63.1;*/
 194     /*par[0] = 63.1;
 195     par[1] = .25;
 196     par[2] = .275;*/
 197     gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
 198
 199     //     Freon
 200     par[0] = geometry->GetOuterFreonWidth()/2;
 201     par[1] = geometry->GetFreonThickness()/2;
 202     par[2] = geometry->GetOuterFreonLength()/2;
 203     /*par[0] = 20.15;
 204     par[1] = .5;                   //Original Settings
 205     par[2] = 65.5;*/
 206     /*par[0] = 65.5;
 207     par[1] = .5;
 208     par[2] = 20.15;*/
 209     gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
 210
 211     par[0] = geometry->GetInnerFreonWidth()/2;
 212     par[1] = geometry->GetFreonThickness()/2;
 213     par[2] = geometry->GetInnerFreonLength()/2;
 214     gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
 215
 216     //     Methane
 217     par[0] = 64.8;
 218     par[1] = geometry->GetGapThickness()/2;
 219     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
 220     par[2] = 64.8;
 221     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
 222
 223     //     Methane gap
 224     par[0] = 64.8;
 225     par[1] = geometry->GetProximityGapThickness()/2;
 226     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
 227     par[2] = 64.8;
 228     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
 229
 230     //     CsI photocathode
 231     par[0] = 64.8;
 232     par[1] = .25;
 233     par[2] = 64.8;
 234     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
 235
 236     //     Anode grid
 237     par[0] = 0.;
 238     par[1] = .001;
 239     par[2] = 20.;
 240     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
 241
 242     // --- Places the detectors defined with GSVOLU
 243     //     Place material inside RICH
 244     gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
 245
 246     gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 -.05 - .376 -.025, 0., 0, "ONLY");
 247     gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2  - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 -.05 - .188, 0., 0, "ONLY");
 248     gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .025, 0., 0, "ONLY");
 249     gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
 250
 251     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
 252
 253     Int_t nspacers = (Int_t)(TMath::Abs(geometry->GetInnerFreonLength()/14.4));
 254     //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers);
 255
 256     //for (i = 1; i <= 9; ++i) {
 257       //zs = (5 - i) * 14.4;                       //Original settings
 258     for (i = 0; i <= nspacers; i++) {
 259         zs = (TMath::Abs(nspacers/2) - i) * 14.4;
 260         gMC->Gspos("SPAC", i, "FRE1", 6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings
 261         //gMC->Gspos("SPAC", i, "FRE1", zs, 0., 6.7, idrotm[1019], "ONLY");
 262     }
 263     //for (i = 10; i <= 18; ++i) {
 264       //zs = (14 - i) * 14.4;                       //Original settings
 265     for (i = nspacers; i < nspacers*2; ++i) {
 266         zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
 267         gMC->Gspos("SPAC", i, "FRE1", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
 268         //gMC->Gspos("SPAC", i, "FRE1", zs, 0., -6.7, idrotm[1019], "ONLY");
 269     }
 270
 271     //for (i = 1; i <= 9; ++i) {
 272       //zs = (5 - i) * 14.4;                       //Original settings
 273       for (i = 0; i <= nspacers; i++) {
 274         zs = (TMath::Abs(nspacers/2) - i) * 14.4;
 275         gMC->Gspos("SPAC", i, "FRE2", 6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings
 276         //gMC->Gspos("SPAC", i, "FRE2", zs, 0., 6.7, idrotm[1019], "ONLY");
 277     }
 278     //for (i = 10; i <= 18; ++i) {
 279       //zs = (5 - i) * 14.4;                       //Original settings
 280       for (i = nspacers; i < nspacers*2; ++i) {
 281         zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
 282         gMC->Gspos("SPAC", i, "FRE2", -6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings
 283         //gMC->Gspos("SPAC", i, "FRE2", zs, 0., -6.7, idrotm[1019], "ONLY");
 284     }
 285
 286     /*gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
 287     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
 288     gMC->Gspos("OQF1", 1, "SRIC", 31.3, -4.724, 41.3, 0, "ONLY");
 289     gMC->Gspos("OQF2", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
 290     gMC->Gspos("OQF1", 3, "SRIC", -31.3, -4.724, -41.3, 0, "ONLY");
 291     gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY");           //Original settings
 292     gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY");            //Original settings
 293     gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
 294     gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
 295     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
 296     gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");*/
 297
 298
 299     gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
 300     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
 301     gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
 302     gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");          //Original settings
 303     gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2), 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");       //Original settings (-31.3)
 304     gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY");           //Original settings
 305     gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY");            //Original settings
 306     gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
 307     gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2+ geometry->GetProximityGapThickness()/2, 0., 0, "ONLY");
 308     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
 309     gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + geometry->GetProximityGapThickness() + .25, 0., 0, "ONLY");
 310
 311     //     Place RICH inside ALICE apparatus
 312
 313     AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
 314     AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
 315     AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
 316     AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
 317     AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
 318     AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
 319     AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
 320
 321     gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
 322     gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
 323     gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
 324     gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
 325     gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
 326     gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
 327     gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
 328
 329 }
 330
 331
 332 //___________________________________________
 333 void AliRICHv0::CreateMaterials()
 334 {
 335     //
 336     // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
 337     // ORIGIN    : NICK VAN EIJNDHOVEN
 338     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
 339     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
 340     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
 341     //
 342     Int_t   ISXFLD = gAlice->Field()->Integ();
 343     Float_t SXMGMX = gAlice->Field()->Max();
 344     Int_t i;
 345
 346     /************************************Antonnelo's Values (14-vectors)*****************************************/
 347     /*
 348     Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
 349                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
 350     Float_t rindex_quarz[14] = { 1.528309,1.533333,
 351                                  1.538243,1.544223,1.550568,1.55777,
 352                                  1.565463,1.574765,1.584831,1.597027,
 353                                1.611858,1.6277,1.6472,1.6724 };
 354     Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 355     Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 356     Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 357     Float_t absco_freon[14] = { 179.0987,179.0987,
 358                                 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
 359     //Float_t absco_freon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
 360         //                       1e-5,1e-5,1e-5,1e-5,1e-5 };
 361     Float_t absco_quarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
 362                                 14.177,9.282,4.0925,1.149,.3627,.10857 };
 363     Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
 364                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
 365     Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
 366                               1e-4,1e-4,1e-4,1e-4 };
 367     Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
 368                                   1e6,1e6,1e6 };
 369     Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
 370                               1e-4,1e-4,1e-4,1e-4 };
 371     Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 372     Float_t effic_csi[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
 373                               .17425,.1785,.1836,.1904,.1938,.221 };
 374     Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 375     */
 376
 377
 378     /**********************************End of Antonnelo's Values**********************************/
 379
 380     /**********************************Values from rich_media.f (31-vectors)**********************************/
 381
 382
 383     //Photons energy intervals
 384     Float_t ppckov[26];
 385     for (i=0;i<26;i++)
 386     {
 387         ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
 388         //printf ("Energy intervals: %e\n",ppckov[i]);
 389     }
 390
 391
 392     //Refraction index for quarz
 393     Float_t rindex_quarz[26];
 394     Float_t  e1= 10.666;
 395     Float_t  e2= 18.125;
 396     Float_t  f1= 46.411;
 397     Float_t  f2= 228.71;
 398     for (i=0;i<26;i++)
 399     {
 400         Float_t ene=ppckov[i]*1e9;
 401         Float_t a=f1/(e1*e1 - ene*ene);
 402         Float_t b=f2/(e2*e2 - ene*ene);
 403         rindex_quarz[i] = TMath::Sqrt(1. + a + b );
 404         //printf ("Rindex_quarz: %e\n",rindex_quarz[i]);
 405     }
 406
 407     //Refraction index for opaque quarz, methane and grid
 408     Float_t rindex_quarzo[26];
 409     Float_t rindex_methane[26];
 410     Float_t rindex_gri[26];
 411     for (i=0;i<26;i++)
 412     {
 413         rindex_quarzo[i]=1;
 414         rindex_methane[i]=1.000444;
 415         rindex_gri[i]=1;
 416         //printf ("Rindex_quarzo , etc: %e, %e, %e\n",rindex_quarzo[i], rindex_methane[i], rindex_gri[i]=1);
 417     }
 418
 419     //Absorption index for freon
 420     Float_t absco_freon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987,  179.0987, 179.0987, 179.0987,
 421                                179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
 422                                7.069031, 4.461292, 2.028366, 1.293013, .577267,   .40746,  .334964, 0., 0., 0.};
 423
 424     //Absorption index for quarz
 425     /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
 426                         .906,.907,.907,.907};
 427     Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
 428                        215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
 429     Float_t absco_quarz[31];
 430     for (Int_t i=0;i<31;i++)
 431     {
 432         Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
 433         if (Xlam <= 160) absco_quarz[i] = 0;
 434         if (Xlam > 250) absco_quarz[i] = 1;
 435         else
 436         {
 437             for (Int_t j=0;j<21;j++)
 438             {
 439                 //printf ("Passed\n");
 440                 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
 441                 {
 442                     Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
 443                     Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
 444                     absco_quarz[i] = -5.0/(TMath::Log(Abso));
 445                 }
 446             }
 447         }
 448         printf ("Absco_quarz: %e Absco_freon: %e for energy: %e\n",absco_quarz[i],absco_freon[i],ppckov[i]);
 449     }*/
 450
 451     /*Float_t absco_quarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
 452                                39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
 453                                17.94531, 11.88753, 5.99128,  3.83503,  2.36661,  1.53155, 1.30582, 1.08574, .8779708,
 454                                .675275, 0., 0., 0.};
 455
 456     for (Int_t i=0;i<31;i++)
 457     {
 458         absco_quarz[i] = absco_quarz[i]/10;
 459     }*/
 460
 461     Float_t absco_quarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
 462                                 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
 463                                 .192, .1497, .10857};
 464
 465     //Absorption index for methane
 466     Float_t absco_methane[26];
 467     for (i=0;i<26;i++)
 468     {
 469         absco_methane[i]=AbsoCH4(ppckov[i]*1e9);
 470         //printf("Absco_methane: %e for energy: %e\n", absco_methane[i],ppckov[i]*1e9);
 471     }
 472
 473     //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
 474     Float_t absco_quarzo[26];
 475     Float_t absco_csi[26];
 476     Float_t absco_gri[26];
 477     Float_t effic_all[26];
 478     Float_t effic_gri[26];
 479     for (i=0;i<26;i++)
 480     {
 481         absco_quarzo[i]=1e-5;
 482         absco_csi[i]=1e-4;
 483         absco_gri[i]=1e-4;
 484         effic_all[i]=1;
 485         effic_gri[i]=1;
 486         //printf ("All must be 1: %e,  %e,  %e,  %e,  %e\n",absco_quarzo[i],absco_csi[i],absco_gri[i],effic_all[i],effic_gri[i]);
 487     }
 488
 489     //Efficiency for csi
 490
 491     Float_t effic_csi[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
 492                              0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
 493                              0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
 494                              0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
 495
 496
 497
 498     //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
 499     //UNPOLARIZED PHOTONS
 500
 501     for (i=0;i<26;i++)
 502     {
 503         effic_csi[i] = effic_csi[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
 504         //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
 505     }
 506
 507     /*******************************************End of rich_media.f***************************************/
 508
 509
 510
 511
 512
 513
 514     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
 515     zmet[2], zqua[2];
 516     Int_t nlmatfre;
 517     Float_t densquao;
 518     Int_t nlmatmet, nlmatqua;
 519     Float_t wmatquao[2], rindex_freon[26];
 520     Float_t aquao[2], epsil, stmin, zquao[2];
 521     Int_t nlmatquao;
 522     Float_t radlal, densal, tmaxfd, deemax, stemax;
 523     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
 524
 525     Int_t *idtmed = fIdtmed->GetArray()-999;
 526
 527     TGeant3 *geant3 = (TGeant3*) gMC;
 528
 529     // --- Photon energy (GeV)
 530     // --- Refraction indexes
 531     for (i = 0; i < 26; ++i) {
 532         rindex_freon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
 533         //printf ("Rindex_freon: %e \n Effic_csi: %e for energy: %e\n",rindex_freon[i], effic_csi[i], ppckov[i]);
 534     }
 535
 536     // --- Detection efficiencies (quantum efficiency for CsI)
 537     // --- Define parameters for honeycomb.
 538     //     Used carbon of equivalent rad. lenght
 539
 540     ahon    = 12.01;
 541     zhon    = 6.;
 542     denshon = 2.265;
 543     radlhon = 18.8;
 544
 545     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
 546
 547     aqua[0]    = 28.09;
 548     aqua[1]    = 16.;
 549     zqua[0]    = 14.;
 550     zqua[1]    = 8.;
 551     densqua    = 2.64;
 552     nlmatqua   = -2;
 553     wmatqua[0] = 1.;
 554     wmatqua[1] = 2.;
 555
 556     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
 557
 558     aquao[0]    = 28.09;
 559     aquao[1]    = 16.;
 560     zquao[0]    = 14.;
 561     zquao[1]    = 8.;
 562     densquao    = 2.64;
 563     nlmatquao   = -2;
 564     wmatquao[0] = 1.;
 565     wmatquao[1] = 2.;
 566
 567     // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
 568
 569     afre[0]    = 12.;
 570     afre[1]    = 19.;
 571     zfre[0]    = 6.;
 572     zfre[1]    = 9.;
 573     densfre    = 1.7;
 574     nlmatfre   = -2;
 575     wmatfre[0] = 6.;
 576     wmatfre[1] = 14.;
 577
 578     // --- Parameters to include in GSMIXT, relative to methane (CH4)
 579
 580     amet[0]    = 12.01;
 581     amet[1]    = 1.;
 582     zmet[0]    = 6.;
 583     zmet[1]    = 1.;
 584     densmet    = 7.17e-4;
 585     nlmatmet   = -2;
 586     wmatmet[0] = 1.;
 587     wmatmet[1] = 4.;
 588
 589     // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
 590
 591     agri    = 63.54;
 592     zgri    = 29.;
 593     densgri = 8.96;
 594     radlgri = 1.43;
 595
 596     // --- Parameters to include in GSMATE related to aluminium sheet
 597
 598     aal    = 26.98;
 599     zal    = 13.;
 600     densal = 2.7;
 601     radlal = 8.9;
 602
 603     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
 604     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
 605     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
 606     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
 607     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
 608     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
 609     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
 610     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
 611     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
 612     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
 613
 614     tmaxfd = -10.;
 615     stemax = -.1;
 616     deemax = -.2;
 617     epsil  = .001;
 618     stmin  = -.001;
 619
 620     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 621     AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 622     AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 623     AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 624     AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 625     AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
 626     AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 627     AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 628     AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
 629     AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 630
 631
 632     geant3->Gsckov(idtmed[1000], 26, ppckov, absco_methane, effic_all, rindex_methane);
 633     geant3->Gsckov(idtmed[1001], 26, ppckov, absco_methane, effic_all, rindex_methane);
 634     geant3->Gsckov(idtmed[1002], 26, ppckov, absco_quarz, effic_all,rindex_quarz);
 635     geant3->Gsckov(idtmed[1003], 26, ppckov, absco_freon, effic_all,rindex_freon);
 636     geant3->Gsckov(idtmed[1004], 26, ppckov, absco_methane, effic_all, rindex_methane);
 637     geant3->Gsckov(idtmed[1005], 26, ppckov, absco_csi, effic_csi, rindex_methane);
 638     geant3->Gsckov(idtmed[1006], 26, ppckov, absco_gri, effic_gri, rindex_gri);
 639     geant3->Gsckov(idtmed[1007], 26, ppckov, absco_quarzo, effic_all, rindex_quarzo);
 640     geant3->Gsckov(idtmed[1008], 26, ppckov, absco_methane, effic_all, rindex_methane);
 641     geant3->Gsckov(idtmed[1009], 26, ppckov, absco_gri, effic_gri, rindex_gri);
 642 }
 643
 644 //___________________________________________
 645
 646 Float_t AliRICHv0::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 647 {
 648
 649     //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
 650
 651     Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
 652                       6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
 653                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
 654
 655
 656     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
 657                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
 658                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
 659                         1.72,1.53};
 660
 661     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
 662                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
 663                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
 664                         1.714,1.498};
 665     Float_t xe=ene;
 666     Int_t  j=Int_t(xe*10)-49;
 667     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
 668     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
 669
 670     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
 671     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 672
 673     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
 674     Float_t tanin=sinin/pdoti;
 675
 676     Float_t c1=cn*cn-ck*ck-sinin*sinin;
 677     Float_t c2=4*cn*cn*ck*ck;
 678     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
 679     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
 680
 681     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
 682     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
 683
 684
 685     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
 686     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
 687
 688     Float_t sigraf=18.;
 689     Float_t lamb=1240/ene;
 690     Float_t fresn;
 691
 692     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
 693
 694     if(pola)
 695     {
 696         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
 697         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
 698     }
 699     else
 700         fresn=0.5*(rp+rs);
 701
 702     fresn = fresn*rO;
 703     return(fresn);
 704 }
 705
 706 //__________________________________________
 707
 708 Float_t AliRICHv0::AbsoCH4(Float_t x)
 709 {
 710
 711     //LOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
 712     Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
 713     //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
 714     Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
 715     const Float_t losch=2.686763E19;                                      // LOSCHMIDT NUMBER IN CM-3
 716     const Float_t igas1=100, igas2=0, oxy=10., wat=5., pre=750.,tem=283.;
 717     Float_t pn=pre/760.;
 718     Float_t tn=tem/273.16;
 719
 720
 721 // ------- METHANE CROSS SECTION -----------------
 722 // ASTROPH. J. 214, L47 (1978)
 723
 724     Float_t sm=0;
 725     if (x<7.75)
 726         sm=.06e-22;
 727
 728     if(x>=7.75 && x<=8.1)
 729     {
 730         Float_t c0=-1.655279e-1;
 731         Float_t c1=6.307392e-2;
 732         Float_t c2=-8.011441e-3;
 733         Float_t c3=3.392126e-4;
 734         sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
 735     }
 736
 737     if (x> 8.1)
 738     {
 739         Int_t j=0;
 740         while (x<=em[j] && x>=em[j+1])
 741         {
 742             j++;
 743             Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
 744             sm=(sch4[j]+a*(x-em[j]))*1e-22;
 745         }
 746     }
 747
 748     Float_t dm=(igas1/100.)*(1.-((oxy+wat)/1.e6))*losch*pn/tn;
 749     Float_t abslm=1./sm/dm;
 750
 751 //    ------- ISOBUTHANE CROSS SECTION --------------
 752 //     i-C4H10 (ai) abs. length from curves in
 753 //     Lu-McDonald paper for BARI RICH workshop .
 754 //     -----------------------------------------------------------
 755
 756     Float_t ai;
 757     Float_t absli;
 758     if (igas2 != 0)
 759     {
 760         if (x<7.25)
 761             ai=100000000.;
 762
 763         if(x>=7.25 && x<7.375)
 764             ai=24.3;
 765
 766         if(x>=7.375)
 767             ai=.0000000001;
 768
 769         Float_t si = 1./(ai*losch*273.16/293.);                    // ISOB. CRO.SEC.IN CM2
 770         Float_t di=(igas2/100.)*(1.-((oxy+wat)/1.e6))*losch*pn/tn;
 771         absli =1./si/di;
 772     }
 773     else
 774         absli=1.e18;
 775 //    ---------------------------------------------------------
 776 //
 777 //       transmission of O2
 778 //
 779 //       y= path in cm, x=energy in eV
 780 //       so= cross section for UV absorption in cm2
 781 //       do= O2 molecular density in cm-3
 782 //    ---------------------------------------------------------
 783
 784     Float_t abslo;
 785     Float_t so=0;
 786     if(x>=6.0)
 787     {
 788         if(x>=6.0 && x<6.5)
 789         {
 790             so=3.392709e-13 * TMath::Exp(2.864104 *x);
 791             so=so*1e-18;
 792         }
 793
 794         if(x>=6.5 && x<7.0)
 795         {
 796             so=2.910039e-34 * TMath::Exp(10.3337*x);
 797             so=so*1e-18;
 798         }
 799
 800
 801         if (x>=7.0)
 802         {
 803             Float_t a0=-73770.76;
 804             Float_t a1=46190.69;
 805             Float_t a2=-11475.44;
 806             Float_t a3=1412.611;
 807             Float_t a4=-86.07027;
 808             Float_t a5=2.074234;
 809             so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
 810             so=so*1e-18;
 811         }
 812
 813         Float_t dox=(oxy/1e6)*losch*pn/tn;
 814         abslo=1./so/dox;
 815     }
 816     else
 817         abslo=1.e18;
 818 //     ---------------------------------------------------------
 819 //
 820 //       transmission of H2O
 821 //
 822 //       y= path in cm, x=energy in eV
 823 //       sw= cross section for UV absorption in cm2
 824 //       dw= H2O molecular density in cm-3
 825 //     ---------------------------------------------------------
 826
 827     Float_t abslw;
 828
 829     Float_t b0=29231.65;
 830     Float_t b1=-15807.74;
 831     Float_t b2=3192.926;
 832     Float_t b3=-285.4809;
 833     Float_t b4=9.533944;
 834
 835     if(x>6.75)
 836     {
 837         Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
 838         sw=sw*1e-18;
 839         Float_t dw=(wat/1e6)*losch*pn/tn;
 840         abslw=1./sw/dw;
 841     }
 842     else
 843         abslw=1.e18;
 844
 845 //    ---------------------------------------------------------
 846
 847     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
 848     return (alength);
 849 }
 850
 851
 852
 853
 854 //___________________________________________
 855
 856 void AliRICHv0::Init()
 857 {
 858     printf("\n\n\n Start Init for version 0 - CPC chamber type \n\n\n");
 859
 860     //
 861     // Initialize Tracking Chambers
 862     //
 863     for (Int_t i=1; i<7; i++) {
 864         //printf ("i:%d",i);
 865         ( (AliRICHChamber*) (*fChambers)[i])->Init();
 866     }
 867
 868     //
 869     // Set the chamber (sensitive region) GEANT identifier
 870
 871     ((AliRICHChamber*)(*fChambers)[0])->SetGid(1);
 872     ((AliRICHChamber*)(*fChambers)[1])->SetGid(2);
 873     ((AliRICHChamber*)(*fChambers)[2])->SetGid(3);
 874     ((AliRICHChamber*)(*fChambers)[3])->SetGid(4);
 875     ((AliRICHChamber*)(*fChambers)[4])->SetGid(5);
 876     ((AliRICHChamber*)(*fChambers)[5])->SetGid(6);
 877     ((AliRICHChamber*)(*fChambers)[6])->SetGid(7);
 878
 879     Float_t pos1[3]={0,471.8999,165.2599};
 880     Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90));
 881
 882     Float_t pos2[3]={171,470,0};
 883     Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90,-20,90,70,0,0));
 884
 885     Float_t pos3[3]={0,500,0};
 886     Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90,0,90,90,0,0));
 887
 888     Float_t pos4[3]={-171,470,0};
 889     Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2], new TRotMatrix("rot996","rot996",90,20,90,110,0,0));
 890
 891     Float_t pos5[3]={161.3999,443.3999,-165.3};
 892     Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70));
 893
 894     Float_t pos6[3]={0., 471.9, -165.3,};
 895     Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90));
 896
 897     Float_t pos7[3]={-161.399,443.3999,-165.3};
 898     Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110));
 899
 900     printf("\n\n\n Finished Init for version 0 - CPC chamber type\n\n\n");
 901 }
 902
 903 //___________________________________________
 904 void AliRICHv0::StepManager()
 905 {
 906     Int_t          copy, id;
 907     static Int_t   idvol;
 908     static Int_t   vol[2];
 909     Int_t          ipart;
 910     static Float_t hits[17];
 911     static Float_t Ckov_data[17];
 912     TLorentzVector Position;
 913     TLorentzVector Momentum;
 914     Float_t        pos[3];
 915     Float_t        mom[4];
 916     Float_t        Localpos[3];
 917     Float_t        Localmom[4];
 918     Float_t        Localtheta,Localphi;
 919     Float_t        theta,phi;
 920     Float_t        destep, step;
 921     Float_t        ranf[2];
 922     static Float_t eloss, xhit, yhit, tlength;
 923     const  Float_t big=1.e10;
 924
 925     TClonesArray &lhits = *fHits;
 926     TGeant3 *geant3 = (TGeant3*) gMC;
 927     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
 928
 929  //if (current->Energy()>1)
 930    //{
 931
 932     // Only gas gap inside chamber
 933     // Tag chambers and record hits when track enters
 934
 935     idvol=-1;
 936     id=gMC->CurrentVolID(copy);
 937     Float_t cherenkov_loss=0;
 938     //gAlice->KeepTrack(gAlice->CurrentTrack());
 939
 940     gMC->TrackPosition(Position);
 941     pos[0]=Position(0);
 942     pos[1]=Position(1);
 943     pos[2]=Position(2);
 944     Ckov_data[1] = pos[0];                 // X-position for hit
 945     Ckov_data[2] = pos[1];                 // Y-position for hit
 946     Ckov_data[3] = pos[2];                 // Z-position for hit
 947     //Ckov_data[11] = gAlice->CurrentTrack();
 948
 949
 950     /********************Store production parameters for Cerenkov photons************************/
 951 //is it a Cerenkov photon?
 952     if (gMC->TrackPid() == 50000050) {
 953
 954       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
 955
 956
 957         Float_t Ckov_energy = current->Energy();
 958         //energy interval for tracking
 959         if  (Ckov_energy > 5.6e-09 && Ckov_energy < 7.8e-09 )
 960         //if (Ckov_energy > 0)
 961         {
 962             if (gMC->IsTrackEntering()){                                     //is track entering?
 963                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 964                   {                                                          //is it in freo?
 965                     if (geant3->Gctrak()->nstep<1){                          //is it the first step?
 966                         Int_t mother = current->GetFirstMother();
 967
 968                         //printf("Second Mother:%d\n",current->GetSecondMother());
 969
 970                         Ckov_data[10] = mother;
 971                         Ckov_data[11] = gAlice->CurrentTrack();
 972                         Ckov_data[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
 973                         fCkov_number++;
 974                         fFreon_prod=1;
 975                         //printf("Index: %d\n",fCkov_number);
 976                     }    //first step question
 977                 }        //freo question
 978
 979                 if (geant3->Gctrak()->nstep<1){                                  //is it first step?
 980                     if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
 981                     {
 982                         Ckov_data[12] = 2;
 983                     }    //quarz question
 984                 }        //first step question
 985
 986                 //printf("Before %d\n",fFreon_prod);
 987             }   //track entering question
 988
 989             if (Ckov_data[12] == 1)                                        //was it produced in Freon?
 990                 //if (fFreon_prod == 1)
 991             {
 992                 if (gMC->IsTrackEntering()){                                     //is track entering?
 993                     //printf("Got in");
 994                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
 995                     {
 996                         //printf("Got in\n");
 997                         gMC->TrackMomentum(Momentum);
 998                         mom[0]=Momentum(0);
 999                         mom[1]=Momentum(1);
1000                         mom[2]=Momentum(2);
1001                         mom[3]=Momentum(3);
1002                                        // Z-position for hit
1003
1004
1005                         /**************** Photons lost in second grid have to be calculated by hand************/
1006
1007                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1008                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
1009                         gMC->Rndm(ranf, 1);
1010                         //printf("grid calculation:%f\n",t);
1011                         if (ranf[0] > t) {
1012                             geant3->StopTrack();
1013                             Ckov_data[13] = 5;
1014                             AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1015                             //printf("Lost one in grid\n");
1016                         }
1017                         /**********************************************************************************/
1018                     }    //gap
1019
1020                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
1021                     {
1022                         gMC->TrackMomentum(Momentum);
1023                         mom[0]=Momentum(0);
1024                         mom[1]=Momentum(1);
1025                         mom[2]=Momentum(2);
1026                         mom[3]=Momentum(3);
1027
1028                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
1029                         /***********************Cerenkov phtons (always polarised)*************************/
1030
1031                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1032                         Float_t t = Fresnel(Ckov_energy*1e9,cophi,1);
1033                         gMC->Rndm(ranf, 1);
1034                         if (ranf[0] < t) {
1035                             geant3->StopTrack();
1036                             Ckov_data[13] = 6;
1037                             AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1038                             //printf("Lost by Fresnel\n");
1039                         }
1040                           /**********************************************************************************/
1041                     }
1042                 } //track entering?
1043
1044
1045                 /********************Evaluation of losses************************/
1046                 /******************still in the old fashion**********************/
1047
1048                 Int_t i1 = geant3->Gctrak()->nmec;            //number of physics mechanisms acting on the particle
1049                 for (Int_t i = 0; i < i1; ++i) {
1050                     //        Reflection loss
1051                     if (geant3->Gctrak()->lmec[i] == 106) {        //was it reflected
1052                         Ckov_data[13]=10;
1053                         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1054                             Ckov_data[13]=1;
1055                         if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1056                             Ckov_data[13]=2;
1057                         geant3->StopTrack();
1058                         AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1059                     } //reflection question
1060
1061
1062                     //        Absorption loss
1063                     else if (geant3->Gctrak()->lmec[i] == 101) {              //was it absorbed?
1064                         Ckov_data[13]=20;
1065                         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1066                             Ckov_data[13]=11;
1067                         if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1068                             Ckov_data[13]=12;
1069                         if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
1070                             Ckov_data[13]=13;
1071                         if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
1072                             Ckov_data[13]=13;
1073
1074                         if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
1075                             Ckov_data[13]=15;
1076
1077                         //        CsI inefficiency
1078                         if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
1079                             Ckov_data[13]=16;
1080                         }
1081                         geant3->StopTrack();
1082                         AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1083                         //printf("Added cerenkov %d\n",fCkov_number);
1084                     } //absorption question
1085
1086
1087                     //        Photon goes out of tracking scope
1088                     else if (geant3->Gctrak()->lmec[i] == 30) {                 //is it below energy treshold?
1089                         Ckov_data[13]=21;
1090                         geant3->StopTrack();
1091                         AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1092                     }   // energy treshold question
1093                 }  //number of mechanisms cycle
1094                 /**********************End of evaluation************************/
1095             } //freon production question
1096         } //energy interval question
1097     } //cerenkov photon question
1098
1099     /**************************************End of Production Parameters Storing*********************/
1100
1101
1102     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
1103
1104     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
1105         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
1106         {
1107
1108           if (gMC->Edep() > 0.){
1109                 gMC->TrackPosition(Position);
1110                 gMC->TrackMomentum(Momentum);
1111                 pos[0]=Position(0);
1112                 pos[1]=Position(1);
1113                 pos[2]=Position(2);
1114                 mom[0]=Momentum(0);
1115                 mom[1]=Momentum(1);
1116                 mom[2]=Momentum(2);
1117                 mom[3]=Momentum(3);
1118                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1119                 Double_t rt = TMath::Sqrt(tc);
1120                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1121                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1122                 gMC->Gmtod(pos,Localpos,1);
1123                 gMC->Gmtod(mom,Localmom,2);
1124
1125                 gMC->CurrentVolOffID(2,copy);
1126                 vol[0]=copy;
1127                 idvol=vol[0]-1;
1128
1129                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1130                         //->Sector(Localpos[0], Localpos[2]);
1131                 //printf("Sector:%d\n",sector);
1132
1133                 /*if (gMC->TrackPid() == 50000051){
1134                   fFeedbacks++;
1135                   printf("Feedbacks:%d\n",fFeedbacks);
1136                 }*/
1137
1138                 ((AliRICHChamber*) (*fChambers)[idvol])
1139                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1140                 if(idvol<7) {
1141                     Ckov_data[0] = gMC->TrackPid();        // particle type
1142                     Ckov_data[1] = pos[0];                 // X-position for hit
1143                     Ckov_data[2] = pos[1];                 // Y-position for hit
1144                     Ckov_data[3] = pos[2];                 // Z-position for hit
1145                     Ckov_data[4] = theta;                      // theta angle of incidence
1146                     Ckov_data[5] = phi;                      // phi angle of incidence
1147                     Ckov_data[8] = (Float_t) fNPadHits;      // first padhit
1148                     Ckov_data[9] = -1;                       // last pad hit
1149                     Ckov_data[13] = 4;                       // photon was detected
1150                     Ckov_data[14] = mom[0];
1151                     Ckov_data[15] = mom[1];
1152                     Ckov_data[16] = mom[2];
1153
1154                     destep = gMC->Edep();
1155                     gMC->SetMaxStep(big);
1156                     cherenkov_loss  += destep;
1157                     Ckov_data[7]=cherenkov_loss;
1158
1159                     MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
1160                     if (fNPadHits > (Int_t)Ckov_data[8]) {
1161                         Ckov_data[8]= Ckov_data[8]+1;
1162                         Ckov_data[9]= (Float_t) fNPadHits;
1163                     }
1164                     //if (sector != -1)
1165                       //{
1166                         AddHit(gAlice->CurrentTrack(),vol,Ckov_data);
1167                         AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1168                       //}
1169                 }
1170             }
1171         }
1172     }
1173
1174     /***********************************************End of photon hits*********************************************/
1175
1176
1177     /**********************************************Charged particles treatment*************************************/
1178
1179     //else if (gMC->TrackCharge())
1180     else if (1 == 1)
1181       {
1182 //If MIP
1183         /*if (gMC->IsTrackEntering())
1184           {
1185             hits[13]=20;//is track entering?
1186           }*/
1187         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1188           {
1189             fFreon_prod=1;
1190           }
1191
1192         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
1193 // Get current particle id (ipart), track position (pos)  and momentum (mom)
1194
1195             gMC->CurrentVolOffID(3,copy);
1196             vol[0]=copy;
1197             idvol=vol[0]-1;
1198
1199             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1200                         //->Sector(Localpos[0], Localpos[2]);
1201             //printf("Sector:%d\n",sector);
1202
1203             gMC->TrackPosition(Position);
1204             gMC->TrackMomentum(Momentum);
1205             pos[0]=Position(0);
1206             pos[1]=Position(1);
1207             pos[2]=Position(2);
1208             mom[0]=Momentum(0);
1209             mom[1]=Momentum(1);
1210             mom[2]=Momentum(2);
1211             mom[3]=Momentum(3);
1212             gMC->Gmtod(pos,Localpos,1);
1213             gMC->Gmtod(mom,Localmom,2);
1214
1215             ipart  = gMC->TrackPid();
1216             //
1217             // momentum loss and steplength in last step
1218             destep = gMC->Edep();
1219             step   = gMC->TrackStep();
1220
1221             //
1222             // record hits when track enters ...
1223             if( gMC->IsTrackEntering()) {
1224 //              gMC->SetMaxStep(fMaxStepGas);
1225                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1226                 Double_t rt = TMath::Sqrt(tc);
1227                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1228                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1229
1230
1231                 Double_t Localtc = Localmom[0]*Localmom[0]+Localmom[2]*Localmom[2];
1232                 Double_t Localrt = TMath::Sqrt(Localtc);
1233                 Localtheta   = Float_t(TMath::ATan2(Localrt,Double_t(Localmom[1])))*kRaddeg;
1234                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[2]),Double_t(Localmom[0])))*kRaddeg;
1235
1236                 hits[0] = Float_t(ipart);         // particle type
1237                 hits[1] = Localpos[0];                 // X-position for hit
1238                 hits[2] = Localpos[1];                 // Y-position for hit
1239                 hits[3] = Localpos[2];                 // Z-position for hit
1240                 hits[4] = Localtheta;                  // theta angle of incidence
1241                 hits[5] = Localphi;                    // phi angle of incidence
1242                 hits[8] = (Float_t) fNPadHits;    // first padhit
1243                 hits[9] = -1;                     // last pad hit
1244                 hits[13] = fFreon_prod;           // did id hit the freon?
1245                 hits[14] = mom[0];
1246                 hits[15] = mom[1];
1247                 hits[16] = mom[2];
1248
1249                 tlength = 0;
1250                 eloss   = 0;
1251                 fFreon_prod = 0;
1252
1253                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
1254
1255
1256                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
1257                 xhit    = Localpos[0];
1258                 yhit    = Localpos[2];
1259                 // Only if not trigger chamber
1260                 if(idvol<7) {
1261                     //
1262                     //  Initialize hit position (cursor) in the segmentation model
1263                     ((AliRICHChamber*) (*fChambers)[idvol])
1264                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1265                 }
1266             }
1267
1268             //
1269             // Calculate the charge induced on a pad (disintegration) in case
1270             //
1271             // Mip left chamber ...
1272             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
1273                 gMC->SetMaxStep(big);
1274                 eloss   += destep;
1275                 tlength += step;
1276
1277
1278                 // Only if not trigger chamber
1279                 if(idvol<7) {
1280                     if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
1281                 }
1282
1283                 hits[6]=tlength;
1284                 hits[7]=eloss;
1285                 if (fNPadHits > (Int_t)hits[8]) {
1286                     hits[8]= hits[8]+1;
1287                     hits[9]= (Float_t) fNPadHits;
1288                 }
1289
1290                 //if(sector !=-1)
1291                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
1292                 eloss = 0;
1293                 //
1294                 // Check additional signal generation conditions
1295                 // defined by the segmentation
1296                 // model (boundary crossing conditions)
1297             } else if
1298                 (((AliRICHChamber*) (*fChambers)[idvol])
1299                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
1300             {
1301                 ((AliRICHChamber*) (*fChambers)[idvol])
1302                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1303                 if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
1304                 xhit     = Localpos[0];
1305                 yhit     = Localpos[2];
1306                 eloss    = destep;
1307                 tlength += step ;
1308                 //
1309                 // nothing special  happened, add up energy loss
1310             } else {
1311                 eloss   += destep;
1312                 tlength += step ;
1313             }
1314         }
1315       }
1316     /*************************************************End of MIP treatment**************************************/
1317    //}
1318 }
1319
1320
1321 //___________________________________________
1322 void AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, Response_t res)
1323 {
1324 //
1325 //  Calls the charge disintegration method of the current chamber and adds
1326 //  the simulated cluster to the root treee
1327 //
1328     Int_t clhits[7];
1329     Float_t newclust[6][500];
1330     Int_t nnew;
1331
1332 //
1333 //  Integrated pulse height on chamber
1334
1335     clhits[0]=fNhits+1;
1336
1337     ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
1338     Int_t ic=0;
1339
1340 //
1341 //  Add new clusters
1342     for (Int_t i=0; i<nnew; i++) {
1343         if (Int_t(newclust[3][i]) > 0) {
1344             ic++;
1345 // Cathode plane
1346             clhits[1] = Int_t(newclust[5][i]);
1347 //  Cluster Charge
1348             clhits[2] = Int_t(newclust[0][i]);
1349 //  Pad: ix
1350             clhits[3] = Int_t(newclust[1][i]);
1351 //  Pad: iy
1352             clhits[4] = Int_t(newclust[2][i]);
1353 //  Pad: charge
1354             clhits[5] = Int_t(newclust[3][i]);
1355 //  Pad: chamber sector
1356             clhits[6] = Int_t(newclust[4][i]);
1357
1358             AddPadHit(clhits);
1359         }
1360     }
1361 }