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