RICH/AliRICHv0.cxx

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