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