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.3  1999/09/29 09:24:29  fca
  19 Introduction of the Copyright and cvs Log
  20
  21 */
  22
  23 /////////////////////////////////////////////////////////
  24 //  Manager and hits classes for set:RICH version 0    //
  25 /////////////////////////////////////////////////////////
  26
  27 #include <TTUBE.h>
  28 #include <TNode.h>
  29 #include <TRandom.h>
  30
  31 #include "AliRICHv0.h"
  32 #include "AliRun.h"
  33 #include "AliMC.h"
  34 #include "iostream.h"
  35 #include "AliCallf77.h"
  36 #include "AliConst.h"
  37 #include "TGeant3.h"
  38
  39 ClassImp(AliRICHv0)
  40
  41 //___________________________________________
  42 AliRICHv0::AliRICHv0() : AliRICH()
  43 {
  44     fChambers = 0;
  45 }
  46
  47 //___________________________________________
  48 AliRICHv0::AliRICHv0(const char *name, const char *title)
  49     : AliRICH(name,title)
  50 {
  51
  52     fChambers = new TObjArray(7);
  53     for (Int_t i=0; i<7; i++) {
  54
  55         (*fChambers)[i] = new AliRICHchamber();
  56
  57     }
  58 }
  59
  60
  61 //___________________________________________
  62 void AliRICHv0::CreateGeometry()
  63 {
  64     //
  65     // Create the geometry for RICH version 1
  66     //
  67     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
  68     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
  69     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
  70     //
  71     //Begin_Html
  72     /*
  73       <img src="picts/AliRICHv1.gif">
  74     */
  75     //End_Html
  76     //Begin_Html
  77     /*
  78       <img src="picts/AliRICHv1Tree.gif">
  79     */
  80     //End_Html
  81
  82
  83     Int_t *idtmed = fIdtmed->GetArray()-999;
  84
  85     Int_t i;
  86     Float_t zs;
  87     Int_t idrotm[1099];
  88     Float_t par[3];
  89
  90     // --- Define the RICH detector
  91     //     External aluminium box
  92     par[0] = 71.1;
  93     par[1] = 11.5;
  94     par[2] = 73.15;
  95     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
  96
  97     //     Sensitive part of the whole RICH
  98     par[0] = 64.8;
  99     par[1] = 11.5;
 100     par[2] = 66.55;
 101     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
 102
 103     //     Honeycomb
 104     par[0] = 63.1;
 105     par[1] = .188;
 106     par[2] = 66.55;
 107     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
 108
 109     //     Aluminium sheet
 110     par[0] = 63.1;
 111     par[1] = .025;
 112     par[2] = 66.55;
 113     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
 114
 115     //     Quartz
 116     par[0] = 63.1;
 117     par[1] = .25;
 118     par[2] = 65.5;
 119     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
 120
 121     //     Spacers (cylinders)
 122     par[0] = 0.;
 123     par[1] = .5;
 124     par[2] = .5;
 125     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
 126
 127     //     Opaque quartz
 128     par[0] = 61.95;
 129     par[1] = .2;
 130     par[2] = 66.5;
 131     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
 132
 133     //     Frame of opaque quartz
 134     par[0] = 20.65;
 135     par[1] = .5;
 136     par[2] = 66.5;
 137     gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
 138
 139     //     Little bar of opaque quartz
 140     par[0] = 63.1;
 141     par[1] = .25;
 142     par[2] = .275;
 143     gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
 144
 145     //     Freon
 146     par[0] = 20.15;
 147     par[1] = .5;
 148     par[2] = 65.5;
 149     gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
 150
 151     //     Methane
 152     par[0] = 64.8;
 153     par[1] = 5.;
 154     par[2] = 64.8;
 155     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
 156
 157     //     Methane gap
 158     par[0] = 64.8;
 159     par[1] = .2;
 160     par[2] = 64.8;
 161     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
 162
 163     //     CsI photocathode
 164     par[0] = 64.8;
 165     par[1] = .25;
 166     par[2] = 64.8;
 167     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
 168
 169     //     Anode grid
 170     par[0] = 0.;
 171     par[1] = .0025;
 172     par[2] = 20.;
 173     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
 174
 175     // --- Places the detectors defined with GSVOLU
 176     //     Place material inside RICH
 177     gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
 178
 179     gMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
 180     gMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
 181     gMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
 182     gMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
 183
 184     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
 185
 186     for (i = 1; i <= 9; ++i) {
 187         zs = (5 - i) * 14.4;
 188         gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
 189     }
 190     for (i = 10; i <= 18; ++i) {
 191         zs = (14 - i) * 14.4;
 192         gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
 193     }
 194
 195     gMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
 196     gMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
 197     gMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
 198     gMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
 199     gMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
 200     gMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
 201     gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
 202     gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
 203     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
 204     gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
 205
 206     //     Place RICH inside ALICE apparatus
 207
 208     AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
 209     AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
 210     AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
 211     AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
 212     AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
 213     AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
 214     AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
 215
 216     gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
 217     gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
 218     gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
 219     gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
 220     gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
 221     gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
 222     gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
 223
 224 }
 225
 226
 227 //___________________________________________
 228 void AliRICHv0::CreateMaterials()
 229 {
 230     //
 231     // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
 232     // ORIGIN    : NICK VAN EIJNDHOVEN
 233     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
 234     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
 235     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
 236     //
 237     Int_t   ISXFLD = gAlice->Field()->Integ();
 238     Float_t SXMGMX = gAlice->Field()->Max();
 239
 240     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,
 241                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
 242     Float_t rindex_quarz[14] = { 1.528309,1.533333,
 243                                  1.538243,1.544223,1.550568,1.55777,
 244                                  1.565463,1.574765,1.584831,1.597027,
 245                                1.611858,1.6277,1.6472,1.6724 };
 246     Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 247     Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 248     Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 249     Float_t absco_freon[14] = { 179.0987,179.0987,
 250                                 179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
 251     Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
 252                                 6.324,4.483,1.6,.323,.073,0. };
 253     Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
 254                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
 255     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,
 256                               1e-4,1e-4,1e-4,1e-4 };
 257     Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
 258                                   1e6,1e6,1e6 };
 259     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,
 260                               1e-4,1e-4,1e-4,1e-4 };
 261     Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 262     Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
 263                               .2101,.2554,.293,.376,.3861,.418 };
 264     Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
 265
 266     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
 267     zmet[2], zqua[2];
 268     Int_t nlmatfre;
 269     Float_t densquao;
 270     Int_t nlmatmet, nlmatqua;
 271     Float_t wmatquao[2], rindex_freon[14];
 272     Int_t i;
 273     Float_t aquao[2], epsil, stmin, zquao[2];
 274     Int_t nlmatquao;
 275     Float_t radlal, densal, tmaxfd, deemax, stemax;
 276     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
 277
 278     Int_t *idtmed = fIdtmed->GetArray()-999;
 279
 280     TGeant3 *geant3 = (TGeant3*) gMC;
 281
 282     // --- Photon energy (GeV)
 283     // --- Refraction indexes
 284     for (i = 0; i < 14; ++i) {
 285         rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
 286     }
 287     // need to be changed
 288
 289     // --- Absorbtion lenghts (in cm)
 290     //      DATA ABSCO_QUARZ /
 291     //     &    5 * 1000000.,
 292     //     &    29.85,    7.34,     4.134,    1.273,    0.722,
 293     //     &    0.365, 0.365, 0.365, 0.  /
 294     // need to be changed
 295
 296     // --- Detection efficiencies (quantum efficiency for CsI)
 297     // --- Define parameters for honeycomb.
 298     //     Used carbon of equivalent rad. lenght
 299
 300     ahon    = 12.01;
 301     zhon    = 6.;
 302     denshon = 2.265;
 303     radlhon = 18.8;
 304
 305     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
 306
 307     aqua[0]    = 28.09;
 308     aqua[1]    = 16.;
 309     zqua[0]    = 14.;
 310     zqua[1]    = 8.;
 311     densqua    = 2.64;
 312     nlmatqua   = -2;
 313     wmatqua[0] = 1.;
 314     wmatqua[1] = 2.;
 315
 316     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
 317
 318     aquao[0]    = 28.09;
 319     aquao[1]    = 16.;
 320     zquao[0]    = 14.;
 321     zquao[1]    = 8.;
 322     densquao    = 2.64;
 323     nlmatquao   = -2;
 324     wmatquao[0] = 1.;
 325     wmatquao[1] = 2.;
 326
 327     // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
 328
 329     afre[0]    = 12.;
 330     afre[1]    = 19.;
 331     zfre[0]    = 6.;
 332     zfre[1]    = 9.;
 333     densfre    = 1.7;
 334     nlmatfre   = -2;
 335     wmatfre[0] = 6.;
 336     wmatfre[1] = 14.;
 337
 338     // --- Parameters to include in GSMIXT, relative to methane (CH4)
 339
 340     amet[0]    = 12.01;
 341     amet[1]    = 1.;
 342     zmet[0]    = 6.;
 343     zmet[1]    = 1.;
 344     densmet    = 7.17e-4;
 345     nlmatmet   = -2;
 346     wmatmet[0] = 1.;
 347     wmatmet[1] = 4.;
 348
 349     // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
 350
 351     agri    = 63.54;
 352     zgri    = 29.;
 353     densgri = 8.96;
 354     radlgri = 1.43;
 355
 356     // --- Parameters to include in GSMATE related to aluminium sheet
 357
 358     aal    = 26.98;
 359     zal    = 13.;
 360     densal = 2.7;
 361     radlal = 8.9;
 362
 363     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
 364     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
 365     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
 366     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
 367     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
 368     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
 369     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
 370     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
 371     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
 372     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
 373
 374     tmaxfd = -10.;
 375     stemax = -.1;
 376     deemax = -.2;
 377     epsil  = .001;
 378     stmin  = -.001;
 379
 380     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 381     AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 382     AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 383     AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 384     AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 385     AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
 386     AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 387     AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 388     AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
 389     AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 390
 391
 392     //     Switch on delta-ray production in the methane and freon gaps
 393
 394     gMC->Gstpar(idtmed[1002], "LOSS", 1.);
 395     gMC->Gstpar(idtmed[1003], "LOSS", 1.);
 396     gMC->Gstpar(idtmed[1004], "LOSS", 1.);
 397     gMC->Gstpar(idtmed[1008], "LOSS", 1.);
 398     gMC->Gstpar(idtmed[1005], "LOSS", 1.);
 399     gMC->Gstpar(idtmed[1002], "HADR", 1.);
 400     gMC->Gstpar(idtmed[1003], "HADR", 1.);
 401     gMC->Gstpar(idtmed[1004], "HADR", 1.);
 402     gMC->Gstpar(idtmed[1008], "HADR", 1.);
 403     gMC->Gstpar(idtmed[1005], "HADR", 1.);
 404     gMC->Gstpar(idtmed[1002], "DCAY", 1.);
 405     gMC->Gstpar(idtmed[1003], "DCAY", 1.);
 406     gMC->Gstpar(idtmed[1004], "DCAY", 1.);
 407     gMC->Gstpar(idtmed[1008], "DCAY", 1.);
 408     gMC->Gstpar(idtmed[1005], "DCAY", 1.);
 409     geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
 410     geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
 411     geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
 412     geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
 413     geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
 414     geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
 415     geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
 416     geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
 417     geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
 418     geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
 419 }
 420
 421 //___________________________________________
 422
 423 void AliRICHv0::Init()
 424 {
 425     printf("\n\n\n Start Init for version 0 - CPC chamber type \n\n\n");
 426
 427     //
 428     // Initialize Tracking Chambers
 429     //
 430     for (Int_t i=0; i<7; i++) {
 431         ( (AliRICHchamber*) (*fChambers)[i])->Init();
 432     }
 433
 434     //
 435     // Set the chamber (sensitive region) GEANT identifier
 436
 437     ((AliRICHchamber*)(*fChambers)[0])->SetGid(1);
 438     ((AliRICHchamber*)(*fChambers)[1])->SetGid(2);
 439     ((AliRICHchamber*)(*fChambers)[2])->SetGid(3);
 440     ((AliRICHchamber*)(*fChambers)[3])->SetGid(4);
 441     ((AliRICHchamber*)(*fChambers)[4])->SetGid(5);
 442     ((AliRICHchamber*)(*fChambers)[5])->SetGid(6);
 443     ((AliRICHchamber*)(*fChambers)[6])->SetGid(7);
 444
 445     printf("\n\n\n Finished Init for version 0 - CPC chamber type\n\n\n");
 446 }
 447
 448 //___________________________________________
 449 void AliRICHv0::StepManager()
 450 {
 451     Int_t          copy, id;
 452     static Int_t   idvol;
 453     static Int_t   vol[2];
 454     Int_t          ipart;
 455     static Float_t hits[10];
 456     TLorentzVector Position;
 457     TLorentzVector Momentum;
 458     Float_t        pos[3];
 459     Float_t        mom[4];
 460     Float_t        Localpos[3];
 461     Float_t        Localmom[4];
 462     Float_t        Localtheta,Localphi;
 463     Float_t        theta,phi;
 464     Float_t        destep, step;
 465     static Float_t eloss, xhit, yhit, tlength;
 466     const  Float_t big=1.e10;
 467
 468     TClonesArray &lhits = *fHits;
 469     TClonesArray &lcerenkovs = *fCerenkovs;
 470
 471     // Only gas gap inside chamber
 472     // Tag chambers and record hits when track enters
 473
 474     idvol=-1;
 475     id=gMC->CurrentVolID(copy);
 476     Float_t cherenkov_loss=0.00001;
 477
 478     //        Treat photons produced in Freon and Quartz
 479     if (gMC->TrackPid() == 50000050 ) {
 480         if (gMC->IsTrackEntering()){
 481             if (gMC->VolId("FREO")==gMC->CurrentVolID(copy) || gMC->VolId("QUAR")==gMC->CurrentVolID(copy)){
 482                 //printf("GOT ONE! Type:%d    \n",gMC->TrackPid());
 483             }
 484         }
 485     }
 486
 487
 488     if (gMC->TrackPid() == 50000050 ) {
 489         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
 490           {
 491             if (gMC->Edep() > 0.){
 492               gMC->TrackPosition(Position);
 493               gMC->TrackMomentum(Momentum);
 494               pos[0]=Position(0);
 495               pos[1]=Position(1);
 496               pos[2]=Position(2);
 497               mom[0]=Momentum(0);
 498               mom[1]=Momentum(1);
 499               mom[2]=Momentum(2);
 500               mom[3]=Momentum(3);
 501               Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 502               Double_t rt = TMath::Sqrt(tc);
 503               theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 504               phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 505               gMC->Gmtod(pos,Localpos,1);
 506               gMC->Gmtod(mom,Localmom,2);
 507
 508               gMC->CurrentVolOffID(2,copy);
 509               vol[0]=copy;
 510               idvol=vol[0]-1;
 511
 512               ((AliRICHchamber*) (*fChambers)[idvol])
 513                 ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 514               if(idvol<7) {
 515                 hits[0] = 50000050;               // particle type
 516                 hits[1] = pos[0];                 // X-position for hit
 517                 hits[2] = pos[1];                 // Y-position for hit
 518                 hits[3] = pos[2];                 // Z-position for hit
 519                 hits[4] = theta;                      // theta angle of incidence
 520                 hits[5] = phi;                      // phi angle of incidence
 521                 hits[8] = (Float_t) fNclusters;   // first padhit
 522                 hits[9] = -1;                     // last pad hit
 523
 524                 MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
 525                 if (fNclusters > (Int_t)hits[8]) {
 526                   hits[8]= hits[8]+1;
 527                   hits[9]= (Float_t) fNclusters;
 528                 }
 529
 530                 AddHit(gAlice->CurrentTrack(),vol,hits);
 531                 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,gAlice->CurrentTrack(),vol,hits);
 532
 533               }
 534             }
 535           }
 536
 537
 538 //
 539 // treat charged particles
 540     } else if (gMC->TrackCharge())
 541     {
 542
 543 //If MIP
 544         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
 545 // Get current particle id (ipart), track position (pos)  and momentum (mom)
 546
 547           gMC->CurrentVolOffID(3,copy);
 548           vol[0]=copy;
 549           idvol=vol[0]-1;
 550
 551             gMC->TrackPosition(Position);
 552             gMC->TrackMomentum(Momentum);
 553             pos[0]=Position(0);
 554             pos[1]=Position(1);
 555             pos[2]=Position(2);
 556             mom[0]=Momentum(0);
 557             mom[1]=Momentum(1);
 558             mom[2]=Momentum(2);
 559             mom[3]=Momentum(3);
 560             gMC->Gmtod(pos,Localpos,1);
 561             gMC->Gmtod(mom,Localmom,2);
 562
 563             ipart  = gMC->TrackPid();
 564             //
 565             // momentum loss and steplength in last step
 566             destep = gMC->Edep();
 567             step   = gMC->TrackStep();
 568
 569             //
 570             // record hits when track enters ...
 571             if( gMC->IsTrackEntering()) {
 572                 gMC->SetMaxStep(fMaxStepGas);
 573                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 574                 Double_t rt = TMath::Sqrt(tc);
 575                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 576                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 577
 578                 Localtheta   = Float_t(TMath::ATan2(rt,Double_t(Localmom[2])))*kRaddeg;
 579                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[1]),Double_t(Localmom[0])))*kRaddeg;
 580
 581                 hits[0] = Float_t(ipart);         // particle type
 582                 hits[1] = pos[0];                 // X-position for hit
 583                 hits[2] = pos[1];                 // Y-position for hit
 584                 hits[3] = pos[2];                 // Z-position for hit
 585                 hits[4] = theta;                  // theta angle of incidence
 586                 hits[5] = phi;                    // phi angle of incidence
 587                 hits[8] = (Float_t) fNclusters;   // first padhit
 588                 hits[9] = -1;                     // last pad hit
 589                 // phi angle of incidence
 590                 tlength = 0;
 591                 eloss   = 0;
 592
 593                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
 594
 595                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
 596                 xhit    = Localpos[0];
 597                 yhit    = Localpos[2];
 598                 // Only if not trigger chamber
 599                 if(idvol<7) {
 600                     //
 601                     //  Initialize hit position (cursor) in the segmentation model
 602                     ((AliRICHchamber*) (*fChambers)[idvol])
 603                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 604                 }
 605             }
 606
 607             //
 608             // Calculate the charge induced on a pad (disintegration) in case
 609             //
 610             // Mip left chamber ...
 611             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 612                 gMC->SetMaxStep(big);
 613                 eloss   += destep;
 614                 tlength += step;
 615
 616
 617
 618                 // Only if not trigger chamber
 619                 if(idvol<7) {
 620                     if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
 621                 }
 622
 623                 hits[6]=tlength;
 624                 hits[7]=eloss;
 625                 if (fNclusters > (Int_t)hits[8]) {
 626                     hits[8]= hits[8]+1;
 627                     hits[9]= (Float_t) fNclusters;
 628                 }
 629
 630                 new(lhits[fNhits++])
 631                     AliRICHhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 632                 eloss = 0;
 633                 //
 634                 // Check additional signal generation conditions
 635                 // defined by the segmentation
 636                 // model (boundary crossing conditions)
 637             } else if
 638                 (((AliRICHchamber*) (*fChambers)[idvol])
 639                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
 640             {
 641                 ((AliRICHchamber*) (*fChambers)[idvol])
 642                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
 643                 if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
 644                 xhit     = Localpos[0];
 645                 yhit     = Localpos[2];
 646                 eloss    = destep;
 647                 tlength += step ;
 648                 //
 649                 // nothing special  happened, add up energy loss
 650             } else {
 651                 eloss   += destep;
 652                 tlength += step ;
 653             }
 654         }
 655     }
 656 }
 657
 658
 659 //___________________________________________
 660 void AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, Response_t res)
 661 {
 662 //
 663 //  Calls the charge disintegration method of the current chamber and adds
 664 //  the simulated cluster to the root treee
 665 //
 666     Int_t clhits[7];
 667     Float_t newclust[6][500];
 668     Int_t nnew;
 669
 670 //
 671 //  Integrated pulse height on chamber
 672
 673     clhits[0]=fNhits+1;
 674
 675     ((AliRICHchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
 676     Int_t ic=0;
 677
 678 //
 679 //  Add new clusters
 680     for (Int_t i=0; i<nnew; i++) {
 681         if (Int_t(newclust[3][i]) > 0) {
 682             ic++;
 683 // Cathode plane
 684             clhits[1] = Int_t(newclust[5][i]);
 685 //  Cluster Charge
 686             clhits[2] = Int_t(newclust[0][i]);
 687 //  Pad: ix
 688             clhits[3] = Int_t(newclust[1][i]);
 689 //  Pad: iy
 690             clhits[4] = Int_t(newclust[2][i]);
 691 //  Pad: charge
 692             clhits[5] = Int_t(newclust[3][i]);
 693 //  Pad: chamber sector
 694             clhits[6] = Int_t(newclust[4][i]);
 695
 696             AddCluster(clhits);
 697         }
 698     }
 699 }
 700
 701 ClassImp(AliRICHchamber)
 702
 703     AliRICHchamber::AliRICHchamber()
 704 {
 705     fSegmentation = new TObjArray(2);
 706     fResponse= new TObjArray(2);
 707     fnsec=1;
 708     frMin=0.1;
 709     frMax=140;
 710 }
 711
 712 //
 713 //  Get reference to response model
 714 AliRICHresponse* AliRICHchamber::GetResponseModel(Response_t res)
 715 {
 716     if (res==mip) {
 717         return (AliRICHresponse*) (*fResponse)[0];
 718     } else if (res==cerenkov) {
 719         return (AliRICHresponse*) (*fResponse)[1];
 720     }
 721     return (AliRICHresponse*) (*fResponse)[0];
 722 }
 723
 724 // Configure response model
 725 void   AliRICHchamber::ResponseModel(Response_t res, AliRICHresponse* thisResponse)
 726 {
 727
 728     if (res==mip) {
 729         (*fResponse)[0]=thisResponse;
 730     } else if (res==cerenkov) {
 731         (*fResponse)[1]=thisResponse;
 732     }
 733 }
 734
 735 void AliRICHchamber::Init()
 736 {
 737
 738     ((AliRICHsegmentation *) (*fSegmentation)[0])->Init(this);
 739     if (fnsec==2) {
 740         ((AliRICHsegmentation *) (*fSegmentation)[1])->Init(this);
 741     }
 742 }
 743
 744 void AliRICHchamber::LocaltoGlobal(Float_t pos[3],Float_t Localpos[3])
 745 {
 746
 747     Double_t *fMatrix;
 748     fMatrix =  fChamberMatrix->GetMatrix();
 749     Localpos[0]=pos[0]*fMatrix[0]+pos[1]*fMatrix[3]+pos[2]*fMatrix[6];
 750     Localpos[1]=pos[0]*fMatrix[1]+pos[1]*fMatrix[4]+pos[2]*fMatrix[7];
 751     Localpos[2]=pos[0]*fMatrix[2]+pos[1]*fMatrix[5]+pos[2]*fMatrix[8];
 752     Localpos[0]+=fChamberTrans[0];
 753     Localpos[1]+=fChamberTrans[1];
 754     Localpos[2]+=fChamberTrans[2];
 755 }
 756
 757
 758 void AliRICHchamber::DisIntegration(Float_t eloss, Float_t xhit, Float_t yhit,
 759                                     Int_t& nnew,Float_t newclust[6][500],Response_t res)
 760 {
 761 //
 762 //  Generates pad hits (simulated cluster)
 763 //  using the segmentation and the response model
 764
 765     Float_t dx, dy;
 766     Float_t local[3];
 767     Float_t source[3];
 768     Int_t Nfp=0;
 769     //    Float_t dummy=0;
 770
 771
 772     //
 773     // Width of the integration area
 774     //
 775     dx=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwX();
 776     dy=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwY();
 777     //
 778     // Get pulse height from energy loss
 779     Float_t qtot;
 780
 781     if (res==mip) {
 782         qtot = ((AliRICHresponse*) (*fResponse)[0])->IntPH(eloss);
 783
 784         local[0]=xhit;
 785         //Z position of the wires relative to the RICH mother volume
 786         local[1]=6.026;
 787         local[2]=yhit;
 788         //Generate feedback photons
 789         Nfp  = ((AliRICHresponse*) (*fResponse)[0])->FeedBackPhotons(source,qtot);
 790         //printf("\nFeedbacks (Mip)     :%d",Nfp);
 791     } else if (res==cerenkov) {
 792         qtot = ((AliRICHresponse*) (*fResponse)[1])->IntPH();
 793         local[0]=xhit;
 794         local[1]=6.026;
 795         local[2]=yhit;
 796         Nfp  = ((AliRICHresponse*) (*fResponse)[1])->FeedBackPhotons(source,qtot);
 797         //printf("\nFeedbacks (Cerenkov):%d",Nfp);
 798     }
 799
 800     //
 801     // Loop Over Pads
 802
 803     Float_t qcheck=0, qp;
 804     nnew=0;
 805     for (Int_t i=1; i<=fnsec; i++) {
 806         qcheck=0;
 807         AliRICHsegmentation * segmentation=(AliRICHsegmentation *) (*fSegmentation)[i-1];
 808         for (segmentation->FirstPad(xhit, yhit, dx, dy);
 809              segmentation->MorePads();
 810              segmentation->NextPad())
 811         {
 812             if (res==mip) {
 813                 qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
 814             }
 815             if (res==cerenkov) {
 816                 qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
 817             }
 818
 819             qp= TMath::Abs(qp);
 820
 821             if (qp > 1.e-4) {
 822                 qcheck+=qp;
 823                 //
 824                 // --- store signal information
 825                 newclust[0][nnew]=qtot;
 826                 newclust[1][nnew]=segmentation->Ix();
 827                 newclust[2][nnew]=segmentation->Iy();
 828                 newclust[3][nnew]=qp * qtot;
 829                 newclust[4][nnew]=segmentation->ISector();
 830                 newclust[5][nnew]=(Float_t) i;
 831                 nnew++;
 832
 833
 834             }
 835         } // Pad loop
 836     } // Cathode plane loop
 837 }
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856