RICH/AliRICHv0.cxx

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