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