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