Corrected error
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
CommitLineData
4c039060 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/*
2e5f0f7b 17 $Log$
5cfcdf54 18 Revision 1.36 2000/12/18 17:45:54 jbarbosa
19 Cleaned up PadHits object.
20
176a9917 21 Revision 1.35 2000/12/15 16:49:40 jbarbosa
22 Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
23
6197ac87 24 Revision 1.34 2000/11/10 18:12:12 jbarbosa
25 Bug fix for AliRICHCerenkov (thanks to P. Hristov)
26
3cc8edee 27 Revision 1.33 2000/11/02 10:09:01 jbarbosa
28 Minor bug correction (some pointers were not initialised in the default constructor)
29
9b674b76 30 Revision 1.32 2000/11/01 15:32:55 jbarbosa
31 Updated to handle both reconstruction algorithms.
32
a4622d0f 33 Revision 1.31 2000/10/26 20:18:33 jbarbosa
34 Supports for methane and freon vessels
35
683b273b 36 Revision 1.30 2000/10/24 13:19:12 jbarbosa
37 Geometry updates.
38
e293afae 39 Revision 1.29 2000/10/19 19:39:25 jbarbosa
40 Some more changes to geometry. Further correction of digitisation "per part. type"
41
53a60579 42 Revision 1.28 2000/10/17 20:50:57 jbarbosa
43 Inversed digtise by particle type (now, only the selected particle type is not digitsed).
44 Corrected several geometry minor bugs.
45 Added new parameter (opaque quartz thickness).
46
3587e146 47 Revision 1.27 2000/10/11 10:33:55 jbarbosa
48 Corrected bug introduced by earlier revisions (CerenkovData array cannot be reset to zero on wach call of StepManager)
49
8cda28a5 50 Revision 1.26 2000/10/03 21:44:08 morsch
51 Use AliSegmentation and AliHit abstract base classes.
52
a2f7eaf6 53 Revision 1.25 2000/10/02 21:28:12 fca
54 Removal of useless dependecies via forward declarations
55
94de3818 56 Revision 1.24 2000/10/02 15:43:17 jbarbosa
57 Fixed forward declarations.
58 Fixed honeycomb density.
59 Fixed cerenkov storing.
60 New electronics.
61
9dda3582 62 Revision 1.23 2000/09/13 10:42:14 hristov
63 Minor corrections for HP, DEC and Sun; strings.h included
64
e813258a 65 Revision 1.22 2000/09/12 18:11:13 fca
66 zero hits area before using
67
e51a3fa9 68 Revision 1.21 2000/07/21 10:21:07 morsch
69 fNrawch = 0; and fNrechits = 0; in the default constructor.
70
db481a6e 71 Revision 1.20 2000/07/10 15:28:39 fca
72 Correction of the inheritance scheme
73
452a64c6 74 Revision 1.19 2000/06/30 16:29:51 dibari
75 Added kDebugLevel variable to control output size on demand
76
33984590 77 Revision 1.18 2000/06/12 15:15:46 jbarbosa
78 Cleaned up version.
79
237c933d 80 Revision 1.17 2000/06/09 14:58:37 jbarbosa
81 New digitisation per particle type
82
d28b5fc2 83 Revision 1.16 2000/04/19 12:55:43 morsch
84 Newly structured and updated version (JB, AM)
85
4c039060 86*/
87
2e5f0f7b 88
ddae0931 89////////////////////////////////////////////////
90// Manager and hits classes for set:RICH //
91////////////////////////////////////////////////
fe4da5cc 92
ddae0931 93#include <TBRIK.h>
94#include <TTUBE.h>
fe4da5cc 95#include <TNode.h>
96#include <TRandom.h>
ddae0931 97#include <TObject.h>
98#include <TVector.h>
99#include <TObjArray.h>
2e5f0f7b 100#include <TArrayF.h>
237c933d 101#include <TFile.h>
102#include <TParticle.h>
94de3818 103#include <TGeometry.h>
104#include <TTree.h>
105
237c933d 106#include <iostream.h>
e813258a 107#include <strings.h>
fe4da5cc 108
fe4da5cc 109#include "AliRICH.h"
a2f7eaf6 110#include "AliSegmentation.h"
237c933d 111#include "AliRICHHit.h"
112#include "AliRICHCerenkov.h"
113#include "AliRICHPadHit.h"
114#include "AliRICHDigit.h"
115#include "AliRICHTransientDigit.h"
116#include "AliRICHRawCluster.h"
a4622d0f 117#include "AliRICHRecHit1D.h"
118#include "AliRICHRecHit3D.h"
237c933d 119#include "AliRICHHitMapA1.h"
2e5f0f7b 120#include "AliRICHClusterFinder.h"
fe4da5cc 121#include "AliRun.h"
ddae0931 122#include "AliMC.h"
94de3818 123#include "AliMagF.h"
452a64c6 124#include "AliConst.h"
125#include "AliPDG.h"
2e5f0f7b 126#include "AliPoints.h"
ddae0931 127#include "AliCallf77.h"
237c933d 128
ddae0931 129
130// Static variables for the pad-hit iterator routines
131static Int_t sMaxIterPad=0;
132static Int_t sCurIterPad=0;
133static TClonesArray *fClusters2;
134static TClonesArray *fHits2;
2e5f0f7b 135static TTree *TrH1;
ddae0931 136
fe4da5cc 137ClassImp(AliRICH)
ddae0931 138
139//___________________________________________
fe4da5cc 140AliRICH::AliRICH()
141{
237c933d 142// Default constructor for RICH manager class
143
ddae0931 144 fIshunt = 0;
145 fHits = 0;
2e5f0f7b 146 fPadHits = 0;
147 fNPadHits = 0;
148 fNcerenkovs = 0;
ddae0931 149 fDchambers = 0;
9b674b76 150 fRecHits1D = 0;
151 fRecHits3D = 0;
152 fRawClusters = 0;
153 fChambers = 0;
ddae0931 154 fCerenkovs = 0;
9dda3582 155 for (Int_t i=0; i<7; i++)
156 {
157 fNdch[i] = 0;
158 fNrawch[i] = 0;
a4622d0f 159 fNrechits1D[i] = 0;
160 fNrechits3D[i] = 0;
9dda3582 161 }
9b674b76 162
163 fFileName = 0;
fe4da5cc 164}
ddae0931 165
166//___________________________________________
fe4da5cc 167AliRICH::AliRICH(const char *name, const char *title)
ddae0931 168 : AliDetector(name,title)
fe4da5cc 169{
ddae0931 170//Begin_Html
171/*
172 <img src="gif/alirich.gif">
173*/
174//End_Html
175
2e5f0f7b 176 fHits = new TClonesArray("AliRICHHit",1000 );
1cedd08a 177 gAlice->AddHitList(fHits);
2e5f0f7b 178 fPadHits = new TClonesArray("AliRICHPadHit",100000);
179 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
180 gAlice->AddHitList(fCerenkovs);
181 //gAlice->AddHitList(fHits);
182 fNPadHits = 0;
183 fNcerenkovs = 0;
184 fIshunt = 0;
ddae0931 185
9dda3582 186 //fNdch = new Int_t[kNCH];
ddae0931 187
237c933d 188 fDchambers = new TObjArray(kNCH);
2e5f0f7b 189
a4622d0f 190 fRecHits1D = new TObjArray(kNCH);
191 fRecHits3D = new TObjArray(kNCH);
ddae0931 192
193 Int_t i;
194
237c933d 195 for (i=0; i<kNCH ;i++) {
2e5f0f7b 196 (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
ddae0931 197 fNdch[i]=0;
198 }
2e5f0f7b 199
9dda3582 200 //fNrawch = new Int_t[kNCH];
ddae0931 201
237c933d 202 fRawClusters = new TObjArray(kNCH);
d28b5fc2 203 //printf("Created fRwClusters with adress:%p",fRawClusters);
2e5f0f7b 204
237c933d 205 for (i=0; i<kNCH ;i++) {
2e5f0f7b 206 (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
207 fNrawch[i]=0;
208 }
209
9dda3582 210 //fNrechits = new Int_t[kNCH];
ddae0931 211
237c933d 212 for (i=0; i<kNCH ;i++) {
a4622d0f 213 (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
214 }
215 for (i=0; i<kNCH ;i++) {
216 (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
2e5f0f7b 217 }
d28b5fc2 218 //printf("Created fRecHits with adress:%p",fRecHits);
2e5f0f7b 219
220
ddae0931 221 SetMarkerColor(kRed);
9b674b76 222
3cc8edee 223 /*fChambers = new TObjArray(kNCH);
9b674b76 224 for (i=0; i<kNCH; i++)
3cc8edee 225 (*fChambers)[i] = new AliRICHChamber();*/
9b674b76 226
227 fFileName = 0;
228
fe4da5cc 229}
230
237c933d 231AliRICH::AliRICH(const AliRICH& RICH)
232{
233// Copy Constructor
234}
235
236
ddae0931 237//___________________________________________
fe4da5cc 238AliRICH::~AliRICH()
239{
237c933d 240
241// Destructor of RICH manager class
242
ddae0931 243 fIshunt = 0;
244 delete fHits;
2e5f0f7b 245 delete fPadHits;
ddae0931 246 delete fCerenkovs;
fe4da5cc 247}
248
ddae0931 249//___________________________________________
fe4da5cc 250void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
251{
237c933d 252
253//
254// Adds a hit to the Hits list
255//
256
ddae0931 257 TClonesArray &lhits = *fHits;
2e5f0f7b 258 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
fe4da5cc 259}
fe4da5cc 260//_____________________________________________________________________________
ddae0931 261void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
262{
237c933d 263
264//
265// Adds a RICH cerenkov hit to the Cerenkov Hits list
266//
267
ddae0931 268 TClonesArray &lcerenkovs = *fCerenkovs;
269 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
2e5f0f7b 270 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
fe4da5cc 271}
ddae0931 272//___________________________________________
2e5f0f7b 273void AliRICH::AddPadHit(Int_t *clhits)
ddae0931 274{
237c933d 275
276//
277// Add a RICH pad hit to the list
278//
279
2e5f0f7b 280 TClonesArray &lPadHits = *fPadHits;
281 new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
33984590 282}
fe4da5cc 283//_____________________________________________________________________________
ddae0931 284void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
285{
237c933d 286
287 //
288 // Add a RICH digit to the list
289 //
ddae0931 290
291 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
2e5f0f7b 292 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
fe4da5cc 293}
294
295//_____________________________________________________________________________
2e5f0f7b 296void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
fe4da5cc 297{
ddae0931 298 //
2e5f0f7b 299 // Add a RICH digit to the list
ddae0931 300 //
2e5f0f7b 301
302 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
303 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
fe4da5cc 304}
305
2e5f0f7b 306//_____________________________________________________________________________
a4622d0f 307void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
308{
309
310 //
311 // Add a RICH reconstructed hit to the list
312 //
313
314 TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
315 new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
316}
317
318//_____________________________________________________________________________
319void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
2e5f0f7b 320{
237c933d 321
322 //
323 // Add a RICH reconstructed hit to the list
324 //
fe4da5cc 325
a4622d0f 326 TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
327 new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
2e5f0f7b 328}
fe4da5cc 329
ddae0931 330//___________________________________________
331void AliRICH::BuildGeometry()
332
fe4da5cc 333{
237c933d 334
335 //
336 // Builds a TNode geometry for event display
337 //
338 TNode *node, *top;
ddae0931 339
340 const int kColorRICH = kGreen;
341 //
237c933d 342 top=gAlice->GetGeometry()->GetNode("alice");
ddae0931 343
344
345 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
346
237c933d 347 top->cd();
ddae0931 348 Float_t pos1[3]={0,471.8999,165.2599};
2e5f0f7b 349 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
350 new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
237c933d 351 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
ddae0931 352
353
237c933d 354 node->SetLineColor(kColorRICH);
355 fNodes->Add(node);
356 top->cd();
ddae0931 357
358 Float_t pos2[3]={171,470,0};
2e5f0f7b 359 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
360 new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
237c933d 361 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
ddae0931 362
363
237c933d 364 node->SetLineColor(kColorRICH);
365 fNodes->Add(node);
366 top->cd();
ddae0931 367 Float_t pos3[3]={0,500,0};
2e5f0f7b 368 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
369 new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
237c933d 370 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
ddae0931 371
372
237c933d 373 node->SetLineColor(kColorRICH);
374 fNodes->Add(node);
375 top->cd();
ddae0931 376 Float_t pos4[3]={-171,470,0};
2e5f0f7b 377 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
378 new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
237c933d 379 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
ddae0931 380
381
237c933d 382 node->SetLineColor(kColorRICH);
383 fNodes->Add(node);
384 top->cd();
ddae0931 385 Float_t pos5[3]={161.3999,443.3999,-165.3};
2e5f0f7b 386 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
387 new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
237c933d 388 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
ddae0931 389
237c933d 390 node->SetLineColor(kColorRICH);
391 fNodes->Add(node);
392 top->cd();
ddae0931 393 Float_t pos6[3]={0., 471.9, -165.3,};
2e5f0f7b 394 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
395 new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
237c933d 396 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
ddae0931 397
398
237c933d 399 node->SetLineColor(kColorRICH);
400 fNodes->Add(node);
401 top->cd();
ddae0931 402 Float_t pos7[3]={-161.399,443.3999,-165.3};
2e5f0f7b 403 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
404 new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
237c933d 405 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
406 node->SetLineColor(kColorRICH);
407 fNodes->Add(node);
ddae0931 408
fe4da5cc 409}
410
452a64c6 411//___________________________________________
412void AliRICH::CreateGeometry()
413{
414 //
415 // Create the geometry for RICH version 1
416 //
417 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
418 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
419 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
420 //
421 //Begin_Html
422 /*
423 <img src="picts/AliRICHv1.gif">
424 */
425 //End_Html
426 //Begin_Html
427 /*
428 <img src="picts/AliRICHv1Tree.gif">
429 */
430 //End_Html
431
432 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
a2f7eaf6 433 AliSegmentation* segmentation;
452a64c6 434 AliRICHGeometry* geometry;
435 AliRICHChamber* iChamber;
436
437 iChamber = &(pRICH->Chamber(0));
438 segmentation=iChamber->GetSegmentationModel(0);
439 geometry=iChamber->GetGeometryModel();
440
441 Float_t distance;
442 distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
443 geometry->SetRadiatorToPads(distance);
444
3587e146 445 //Opaque quartz thickness
e293afae 446 Float_t oqua_thickness = .5;
683b273b 447 //CsI dimensions
448
449 Float_t csi_length = 160*.8 + 2.6;
450 Float_t csi_width = 144*.84 + 2*2.6;
452a64c6 451
452 Int_t *idtmed = fIdtmed->GetArray()-999;
453
454 Int_t i;
455 Float_t zs;
456 Int_t idrotm[1099];
457 Float_t par[3];
458
459 // --- Define the RICH detector
460 // External aluminium box
e293afae 461 par[0] = 68.8;
6197ac87 462 par[1] = 13; //Original Settings
e293afae 463 par[2] = 70.86;
452a64c6 464 /*par[0] = 73.15;
465 par[1] = 11.5;
466 par[2] = 71.1;*/
467 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
468
683b273b 469 // Air
470 par[0] = 66.3;
6197ac87 471 par[1] = 13; //Original Settings
683b273b 472 par[2] = 68.35;
452a64c6 473 /*par[0] = 66.55;
474 par[1] = 11.5;
475 par[2] = 64.8;*/
476 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
477
6197ac87 478 // Air 2 (cutting the lower part of the box)
479
480 par[0] = 1.25;
481 par[1] = 3; //Original Settings
482 par[2] = 70.86;
483 gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
484
485 // Air 3 (cutting the lower part of the box)
486
487 par[0] = 66.3;
488 par[1] = 3; //Original Settings
489 par[2] = 1.2505;
490 gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
491
452a64c6 492 // Honeycomb
e293afae 493 par[0] = 66.3;
452a64c6 494 par[1] = .188; //Original Settings
e293afae 495 par[2] = 68.35;
452a64c6 496 /*par[0] = 66.55;
497 par[1] = .188;
498 par[2] = 63.1;*/
499 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
500
501 // Aluminium sheet
e293afae 502 par[0] = 66.3;
452a64c6 503 par[1] = .025; //Original Settings
e293afae 504 par[2] = 68.35;
452a64c6 505 /*par[0] = 66.5;
506 par[1] = .025;
507 par[2] = 63.1;*/
508 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
509
510 // Quartz
511 par[0] = geometry->GetQuartzWidth()/2;
512 par[1] = geometry->GetQuartzThickness()/2;
513 par[2] = geometry->GetQuartzLength()/2;
514 /*par[0] = 63.1;
515 par[1] = .25; //Original Settings
516 par[2] = 65.5;*/
517 /*par[0] = geometry->GetQuartzWidth()/2;
518 par[1] = geometry->GetQuartzThickness()/2;
519 par[2] = geometry->GetQuartzLength()/2;*/
520 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
521 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
522
523 // Spacers (cylinders)
524 par[0] = 0.;
525 par[1] = .5;
526 par[2] = geometry->GetFreonThickness()/2;
527 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
528
683b273b 529 // Feet (freon slabs supports)
530
531 par[0] = .7;
532 par[1] = .3;
533 par[2] = 1.9;
534 gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
535
452a64c6 536 // Opaque quartz
e293afae 537 par[0] = geometry->GetQuartzWidth()/2;
538 par[1] = .2;
539 par[2] = geometry->GetQuartzLength()/2;
540 /*par[0] = 61.95;
452a64c6 541 par[1] = .2; //Original Settings
e293afae 542 par[2] = 66.5;*/
452a64c6 543 /*par[0] = 66.5;
544 par[1] = .2;
545 par[2] = 61.95;*/
546 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
547
548 // Frame of opaque quartz
53a60579 549 par[0] = geometry->GetOuterFreonWidth()/2;
550 //+ oqua_thickness;
452a64c6 551 par[1] = geometry->GetFreonThickness()/2;
53a60579 552 par[2] = geometry->GetOuterFreonLength()/2;
553 //+ oqua_thickness;
452a64c6 554 /*par[0] = 20.65;
555 par[1] = .5; //Original Settings
556 par[2] = 66.5;*/
557 /*par[0] = 66.5;
558 par[1] = .5;
559 par[2] = 20.65;*/
560 gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
561
e293afae 562 par[0] = geometry->GetInnerFreonWidth()/2;
452a64c6 563 par[1] = geometry->GetFreonThickness()/2;
e293afae 564 par[2] = geometry->GetInnerFreonLength()/2;
452a64c6 565 gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
566
567 // Little bar of opaque quartz
e293afae 568 //par[0] = .275;
569 //par[1] = geometry->GetQuartzThickness()/2;
53a60579 570 //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
e293afae 571 //par[2] = geometry->GetInnerFreonLength()/2;
53a60579 572 //+ oqua_thickness;
452a64c6 573 /*par[0] = .275;
574 par[1] = .25; //Original Settings
575 par[2] = 63.1;*/
576 /*par[0] = 63.1;
577 par[1] = .25;
578 par[2] = .275;*/
e293afae 579 //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
452a64c6 580
581 // Freon
e293afae 582 par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
452a64c6 583 par[1] = geometry->GetFreonThickness()/2;
e293afae 584 par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
452a64c6 585 /*par[0] = 20.15;
586 par[1] = .5; //Original Settings
587 par[2] = 65.5;*/
588 /*par[0] = 65.5;
589 par[1] = .5;
590 par[2] = 20.15;*/
591 gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
592
e293afae 593 par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
452a64c6 594 par[1] = geometry->GetFreonThickness()/2;
e293afae 595 par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
452a64c6 596 gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
597
598 // Methane
e293afae 599 //par[0] = 64.8;
683b273b 600 par[0] = csi_width/2;
452a64c6 601 par[1] = geometry->GetGapThickness()/2;
602 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
e293afae 603 //par[2] = 64.8;
683b273b 604 par[2] = csi_length/2;
452a64c6 605 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
606
607 // Methane gap
e293afae 608 //par[0] = 64.8;
683b273b 609 par[0] = csi_width/2;
452a64c6 610 par[1] = geometry->GetProximityGapThickness()/2;
611 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
e293afae 612 //par[2] = 64.8;
683b273b 613 par[2] = csi_length/2;
452a64c6 614 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
615
616 // CsI photocathode
e293afae 617 //par[0] = 64.8;
683b273b 618 par[0] = csi_width/2;
452a64c6 619 par[1] = .25;
e293afae 620 //par[2] = 64.8;
683b273b 621 par[2] = csi_length/2;
452a64c6 622 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
623
624 // Anode grid
625 par[0] = 0.;
626 par[1] = .001;
627 par[2] = 20.;
628 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
683b273b 629
6197ac87 630 // Wire supports
631 // Bar of metal
632
633 par[0] = csi_width/2;
634 par[1] = 1.05;
635 par[2] = 1.05;
636 gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
637
638 // Ceramic pick up (base)
639
640 par[0] = csi_width/2;
641 par[1] = .25;
642 par[2] = 1.05;
643 gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
644
645 // Ceramic pick up (head)
646
647 par[0] = csi_width/2;
648 par[1] = .1;
649 par[2] = .1;
650 gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
651
683b273b 652 // Aluminium supports for methane and CsI
653 // Short bar
654
655 par[0] = csi_width/2;
656 par[1] = geometry->GetGapThickness()/2 + .25;
657 par[2] = (68.35 - csi_length/2)/2;
658 gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
452a64c6 659
683b273b 660 // Long bar
661
662 par[0] = (66.3 - csi_width/2)/2;
663 par[1] = geometry->GetGapThickness()/2 + .25;
664 par[2] = csi_length/2 + 68.35 - csi_length/2;
665 gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
666
667 // Aluminium supports for freon
668 // Short bar
669
670 par[0] = geometry->GetQuartzWidth()/2;
671 par[1] = .3;
672 par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
673 gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
674
675 // Long bar
676
677 par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
678 par[1] = .3;
679 par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
680 gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
681
6197ac87 682 // PCB backplane
683b273b 683
6197ac87 684 par[0] = csi_width/2;
685 par[1] = .25;
686 par[2] = csi_length/4 -.5025;
687 gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
683b273b 688
6197ac87 689
690 // Backplane supports
691
692 // Aluminium slab
693
694 par[0] = 33.15;
695 par[1] = 2;
696 par[2] = 21.65;
697 gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
698
699 // Big hole
700
701 par[0] = 9.05;
702 par[1] = 2;
703 par[2] = 4.4625;
704 gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
705
706 // Small hole
707
708 par[0] = 5.7;
709 par[1] = 2;
710 par[2] = 4.4625;
711 gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
712
713 // Place holes inside backplane support
714
715 gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
716 gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
717 gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
718 gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
719 gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
720 gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
721 gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
722 gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
723 gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
724 gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
725 gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
726 gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
727 gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
728 gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
729 gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
730 gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
731
732
733
452a64c6 734 // --- Places the detectors defined with GSVOLU
735 // Place material inside RICH
6197ac87 736 gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
737 gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
738 gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
739 gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, -68.35 - 1.25, 0, "ONLY");
740 gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 68.35 + 1.25, 0, "ONLY");
741
683b273b 742
743 gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
744 gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
745 gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
746 gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
747 gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
748 gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
749 gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
750 gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
751 gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
752 gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
753 gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
452a64c6 754 gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
755
683b273b 756 // Supports placing
757
758 // Methane supports
759 gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
760 gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
761 gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
762 gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
763
764 //Freon supports
765
766 Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
767
768 gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
769 gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
770 gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
771 gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
772
452a64c6 773 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
774
e293afae 775 //Int_t nspacers = (Int_t)(TMath::Abs(geometry->GetInnerFreonLength()/14.4));
776 Int_t nspacers = 9;
452a64c6 777 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers);
778
779 //printf("Nspacers: %d", nspacers);
780
781 //for (i = 1; i <= 9; ++i) {
782 //zs = (5 - i) * 14.4; //Original settings
783 for (i = 0; i < nspacers; i++) {
784 zs = (TMath::Abs(nspacers/2) - i) * 14.4;
785 gMC->Gspos("SPAC", i, "FRE1", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
786 //gMC->Gspos("SPAC", i, "FRE1", zs, 0., 6.7, idrotm[1019], "ONLY");
787 }
788 //for (i = 10; i <= 18; ++i) {
789 //zs = (14 - i) * 14.4; //Original settings
790 for (i = nspacers; i < nspacers*2; ++i) {
791 zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
792 gMC->Gspos("SPAC", i, "FRE1", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
793 //gMC->Gspos("SPAC", i, "FRE1", zs, 0., -6.7, idrotm[1019], "ONLY");
794 }
795
796 //for (i = 1; i <= 9; ++i) {
797 //zs = (5 - i) * 14.4; //Original settings
798 for (i = 0; i < nspacers; i++) {
799 zs = (TMath::Abs(nspacers/2) - i) * 14.4;
800 gMC->Gspos("SPAC", i, "FRE2", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
801 //gMC->Gspos("SPAC", i, "FRE2", zs, 0., 6.7, idrotm[1019], "ONLY");
802 }
803 //for (i = 10; i <= 18; ++i) {
804 //zs = (5 - i) * 14.4; //Original settings
805 for (i = nspacers; i < nspacers*2; ++i) {
806 zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
807 gMC->Gspos("SPAC", i, "FRE2", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
808 //gMC->Gspos("SPAC", i, "FRE2", zs, 0., -6.7, idrotm[1019], "ONLY");
809 }
810
811 /*gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
812 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
813 gMC->Gspos("OQF1", 1, "SRIC", 31.3, -4.724, 41.3, 0, "ONLY");
814 gMC->Gspos("OQF2", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
815 gMC->Gspos("OQF1", 3, "SRIC", -31.3, -4.724, -41.3, 0, "ONLY");
816 gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY"); //Original settings
817 gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY"); //Original settings
818 gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
819 gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
820 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
821 gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");*/
822
823
824 gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
825 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
683b273b 826 gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
452a64c6 827 gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings
683b273b 828 gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3)
e293afae 829 //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65)
830 //gMC->Gspos("BARR", 2, "QUAR", geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (21.65)
452a64c6 831 gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
832 gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
833 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
834 gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
6197ac87 835
836 // Wire support placing
837
838 gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
839 gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
840 gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
841
842 // Backplane placing
843
844 gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
845 gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
846 gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
847 gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
848 gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
849 gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
850
851 // PCB placing
852
853 gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 1.2, 0, "ONLY");
854 gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 1.2, 0, "ONLY");
855
856
452a64c6 857
858 //printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2);
859
860 // Place RICH inside ALICE apparatus
861
862 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
863 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
864 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
865 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
866 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
867 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
868 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
869
870 gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
871 gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
872 gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
873 gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
874 gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
875 gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
876 gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
877
878}
879
880
881//___________________________________________
882void AliRICH::CreateMaterials()
883{
884 //
885 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
886 // ORIGIN : NICK VAN EIJNDHOVEN
887 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
888 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
889 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
890 //
891 Int_t isxfld = gAlice->Field()->Integ();
892 Float_t sxmgmx = gAlice->Field()->Max();
893 Int_t i;
894
895 /************************************Antonnelo's Values (14-vectors)*****************************************/
896 /*
897 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,
898 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
899 Float_t rIndexQuarz[14] = { 1.528309,1.533333,
900 1.538243,1.544223,1.550568,1.55777,
901 1.565463,1.574765,1.584831,1.597027,
902 1.611858,1.6277,1.6472,1.6724 };
903 Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
904 Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
905 Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
906 Float_t abscoFreon[14] = { 179.0987,179.0987,
907 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
908 //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
909 // 1e-5,1e-5,1e-5,1e-5,1e-5 };
910 Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
911 14.177,9.282,4.0925,1.149,.3627,.10857 };
912 Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
913 1e-5,1e-5,1e-5,1e-5,1e-5 };
914 Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
915 1e-4,1e-4,1e-4,1e-4 };
916 Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
917 1e6,1e6,1e6 };
918 Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
919 1e-4,1e-4,1e-4,1e-4 };
920 Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
921 Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
922 .17425,.1785,.1836,.1904,.1938,.221 };
923 Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
924 */
925
926
927 /**********************************End of Antonnelo's Values**********************************/
928
929 /**********************************Values from rich_media.f (31-vectors)**********************************/
930
931
932 //Photons energy intervals
933 Float_t ppckov[26];
934 for (i=0;i<26;i++)
935 {
936 ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
937 //printf ("Energy intervals: %e\n",ppckov[i]);
938 }
939
940
941 //Refraction index for quarz
942 Float_t rIndexQuarz[26];
943 Float_t e1= 10.666;
944 Float_t e2= 18.125;
945 Float_t f1= 46.411;
946 Float_t f2= 228.71;
947 for (i=0;i<26;i++)
948 {
949 Float_t ene=ppckov[i]*1e9;
950 Float_t a=f1/(e1*e1 - ene*ene);
951 Float_t b=f2/(e2*e2 - ene*ene);
952 rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
953 //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
954 }
955
956 //Refraction index for opaque quarz, methane and grid
957 Float_t rIndexOpaqueQuarz[26];
958 Float_t rIndexMethane[26];
959 Float_t rIndexGrid[26];
960 for (i=0;i<26;i++)
961 {
962 rIndexOpaqueQuarz[i]=1;
963 rIndexMethane[i]=1.000444;
964 rIndexGrid[i]=1;
965 //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
966 }
967
968 //Absorption index for freon
969 Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
970 179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
971 7.069031, 4.461292, 2.028366, 1.293013, .577267, .40746, .334964, 0., 0., 0.};
972
973 //Absorption index for quarz
974 /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
975 .906,.907,.907,.907};
976 Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
977 215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
978 Float_t abscoQuarz[31];
979 for (Int_t i=0;i<31;i++)
980 {
981 Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
982 if (Xlam <= 160) abscoQuarz[i] = 0;
983 if (Xlam > 250) abscoQuarz[i] = 1;
984 else
985 {
986 for (Int_t j=0;j<21;j++)
987 {
988 //printf ("Passed\n");
989 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
990 {
991 Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
992 Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
993 abscoQuarz[i] = -5.0/(TMath::Log(Abso));
994 }
995 }
996 }
997 printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
998 }*/
999
1000 /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1001 39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1002 17.94531, 11.88753, 5.99128, 3.83503, 2.36661, 1.53155, 1.30582, 1.08574, .8779708,
1003 .675275, 0., 0., 0.};
1004
1005 for (Int_t i=0;i<31;i++)
1006 {
1007 abscoQuarz[i] = abscoQuarz[i]/10;
1008 }*/
1009
1010 Float_t abscoQuarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
1011 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1012 .192, .1497, .10857};
1013
1014 //Absorption index for methane
1015 Float_t abscoMethane[26];
1016 for (i=0;i<26;i++)
1017 {
1018 abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1019 //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1020 }
1021
1022 //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1023 Float_t abscoOpaqueQuarz[26];
1024 Float_t abscoCsI[26];
1025 Float_t abscoGrid[26];
1026 Float_t efficAll[26];
1027 Float_t efficGrid[26];
1028 for (i=0;i<26;i++)
1029 {
1030 abscoOpaqueQuarz[i]=1e-5;
1031 abscoCsI[i]=1e-4;
1032 abscoGrid[i]=1e-4;
1033 efficAll[i]=1;
1034 efficGrid[i]=1;
1035 //printf ("All must be 1: %e, %e, %e, %e, %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1036 }
1037
1038 //Efficiency for csi
1039
1040 Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1041 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1042 0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1043 0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1044
1045
1046
1047 //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1048 //UNPOLARIZED PHOTONS
1049
1050 for (i=0;i<26;i++)
1051 {
1052 efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1053 //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1054 }
1055
1056 /*******************************************End of rich_media.f***************************************/
1057
1058
1059
1060
1061
1062
1063 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1064 zmet[2], zqua[2];
1065 Int_t nlmatfre;
1066 Float_t densquao;
1067 Int_t nlmatmet, nlmatqua;
1068 Float_t wmatquao[2], rIndexFreon[26];
1069 Float_t aquao[2], epsil, stmin, zquao[2];
1070 Int_t nlmatquao;
1071 Float_t radlal, densal, tmaxfd, deemax, stemax;
1072 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1073
1074 Int_t *idtmed = fIdtmed->GetArray()-999;
1075
452a64c6 1076 // --- Photon energy (GeV)
1077 // --- Refraction indexes
1078 for (i = 0; i < 26; ++i) {
1079 rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1080 //rIndexFreon[i] = 1;
1081 //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1082 }
1083
1084 // --- Detection efficiencies (quantum efficiency for CsI)
1085 // --- Define parameters for honeycomb.
1086 // Used carbon of equivalent rad. lenght
1087
1088 ahon = 12.01;
1089 zhon = 6.;
9dda3582 1090 denshon = 0.1;
452a64c6 1091 radlhon = 18.8;
1092
1093 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1094
1095 aqua[0] = 28.09;
1096 aqua[1] = 16.;
1097 zqua[0] = 14.;
1098 zqua[1] = 8.;
1099 densqua = 2.64;
1100 nlmatqua = -2;
1101 wmatqua[0] = 1.;
1102 wmatqua[1] = 2.;
1103
1104 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1105
1106 aquao[0] = 28.09;
1107 aquao[1] = 16.;
1108 zquao[0] = 14.;
1109 zquao[1] = 8.;
1110 densquao = 2.64;
1111 nlmatquao = -2;
1112 wmatquao[0] = 1.;
1113 wmatquao[1] = 2.;
1114
1115 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1116
1117 afre[0] = 12.;
1118 afre[1] = 19.;
1119 zfre[0] = 6.;
1120 zfre[1] = 9.;
1121 densfre = 1.7;
1122 nlmatfre = -2;
1123 wmatfre[0] = 6.;
1124 wmatfre[1] = 14.;
1125
1126 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1127
1128 amet[0] = 12.01;
1129 amet[1] = 1.;
1130 zmet[0] = 6.;
1131 zmet[1] = 1.;
1132 densmet = 7.17e-4;
1133 nlmatmet = -2;
1134 wmatmet[0] = 1.;
1135 wmatmet[1] = 4.;
1136
1137 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1138
1139 agri = 63.54;
1140 zgri = 29.;
1141 densgri = 8.96;
1142 radlgri = 1.43;
1143
1144 // --- Parameters to include in GSMATE related to aluminium sheet
1145
1146 aal = 26.98;
1147 zal = 13.;
1148 densal = 2.7;
1149 radlal = 8.9;
6197ac87 1150
1151 // --- Glass parameters
1152
1153 Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
1154 Float_t zglass[5]={ 6., 14., 8., 5., 11.};
1155 Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
1156 Float_t dglass=1.74;
1157
452a64c6 1158
1159 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1160 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1161 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1162 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1163 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1164 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1165 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1166 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1167 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1168 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
6197ac87 1169 AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1170 AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
452a64c6 1171
1172 tmaxfd = -10.;
1173 stemax = -.1;
1174 deemax = -.2;
1175 epsil = .001;
1176 stmin = -.001;
1177
1178 AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1179 AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1180 AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1181 AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1182 AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1183 AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1184 AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1185 AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1186 AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1187 AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
6197ac87 1188 AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1189 AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
452a64c6 1190
1191
5cfcdf54 1192 gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1193 gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1194 gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1195 gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1196 gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1197 gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1198 gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1199 gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1200 gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1201 gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1202 gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
452a64c6 1203}
1204
1205//___________________________________________
1206
1207Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1208{
1209
1210 //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1211
1212 Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
1213 6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
1214 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1215
1216
1217 Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1218 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1219 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1220 1.72,1.53};
1221
1222 Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1223 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1224 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1225 1.714,1.498};
1226 Float_t xe=ene;
1227 Int_t j=Int_t(xe*10)-49;
1228 Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1229 Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1230
1231 //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1232 //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1233
1234 Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1235 Float_t tanin=sinin/pdoti;
1236
1237 Float_t c1=cn*cn-ck*ck-sinin*sinin;
1238 Float_t c2=4*cn*cn*ck*ck;
1239 Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1240 Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1241
1242 Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1243 Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1244
1245
1246 //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1247 //B.J. STAGG APPLIED OPTICS, 30(1991),4113
1248
1249 Float_t sigraf=18.;
1250 Float_t lamb=1240/ene;
1251 Float_t fresn;
1252
1253 Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1254
1255 if(pola)
1256 {
1257 Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S
1258 fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1259 }
1260 else
1261 fresn=0.5*(rp+rs);
1262
1263 fresn = fresn*rO;
1264 return(fresn);
1265}
1266
1267//__________________________________________
1268Float_t AliRICH::AbsoCH4(Float_t x)
1269{
1270
1271 //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1272 Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22
1273 //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1274 Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1275 const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3
1276 const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1277 Float_t pn=kPressure/760.;
1278 Float_t tn=kTemperature/273.16;
1279
1280
1281// ------- METHANE CROSS SECTION -----------------
1282// ASTROPH. J. 214, L47 (1978)
1283
1284 Float_t sm=0;
1285 if (x<7.75)
1286 sm=.06e-22;
1287
1288 if(x>=7.75 && x<=8.1)
1289 {
1290 Float_t c0=-1.655279e-1;
1291 Float_t c1=6.307392e-2;
1292 Float_t c2=-8.011441e-3;
1293 Float_t c3=3.392126e-4;
1294 sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1295 }
1296
1297 if (x> 8.1)
1298 {
1299 Int_t j=0;
1300 while (x<=em[j] && x>=em[j+1])
1301 {
1302 j++;
1303 Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1304 sm=(sch4[j]+a*(x-em[j]))*1e-22;
1305 }
1306 }
1307
1308 Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1309 Float_t abslm=1./sm/dm;
1310
1311// ------- ISOBUTHANE CROSS SECTION --------------
1312// i-C4H10 (ai) abs. length from curves in
1313// Lu-McDonald paper for BARI RICH workshop .
1314// -----------------------------------------------------------
1315
1316 Float_t ai;
1317 Float_t absli;
1318 if (kIgas2 != 0)
1319 {
1320 if (x<7.25)
1321 ai=100000000.;
1322
1323 if(x>=7.25 && x<7.375)
1324 ai=24.3;
1325
1326 if(x>=7.375)
1327 ai=.0000000001;
1328
1329 Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2
1330 Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1331 absli =1./si/di;
1332 }
1333 else
1334 absli=1.e18;
1335// ---------------------------------------------------------
1336//
1337// transmission of O2
1338//
1339// y= path in cm, x=energy in eV
1340// so= cross section for UV absorption in cm2
1341// do= O2 molecular density in cm-3
1342// ---------------------------------------------------------
1343
1344 Float_t abslo;
1345 Float_t so=0;
1346 if(x>=6.0)
1347 {
1348 if(x>=6.0 && x<6.5)
1349 {
1350 so=3.392709e-13 * TMath::Exp(2.864104 *x);
1351 so=so*1e-18;
1352 }
1353
1354 if(x>=6.5 && x<7.0)
1355 {
1356 so=2.910039e-34 * TMath::Exp(10.3337*x);
1357 so=so*1e-18;
1358 }
1359
1360
1361 if (x>=7.0)
1362 {
1363 Float_t a0=-73770.76;
1364 Float_t a1=46190.69;
1365 Float_t a2=-11475.44;
1366 Float_t a3=1412.611;
1367 Float_t a4=-86.07027;
1368 Float_t a5=2.074234;
1369 so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1370 so=so*1e-18;
1371 }
1372
1373 Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1374 abslo=1./so/dox;
1375 }
1376 else
1377 abslo=1.e18;
1378// ---------------------------------------------------------
1379//
1380// transmission of H2O
1381//
1382// y= path in cm, x=energy in eV
1383// sw= cross section for UV absorption in cm2
1384// dw= H2O molecular density in cm-3
1385// ---------------------------------------------------------
1386
1387 Float_t abslw;
1388
1389 Float_t b0=29231.65;
1390 Float_t b1=-15807.74;
1391 Float_t b2=3192.926;
1392 Float_t b3=-285.4809;
1393 Float_t b4=9.533944;
1394
1395 if(x>6.75)
1396 {
1397 Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1398 sw=sw*1e-18;
1399 Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1400 abslw=1./sw/dw;
1401 }
1402 else
1403 abslw=1.e18;
1404
1405// ---------------------------------------------------------
1406
1407 Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1408 return (alength);
1409}
1410
1411
1412
ddae0931 1413//___________________________________________
1414Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1415{
237c933d 1416
1417// Default value
1418
ddae0931 1419 return 9999;
1420}
fe4da5cc 1421
ddae0931 1422//___________________________________________
1423void AliRICH::MakeBranch(Option_t* option)
fe4da5cc 1424{
237c933d 1425 // Create Tree branches for the RICH.
fe4da5cc 1426
237c933d 1427 const Int_t kBufferSize = 4000;
ddae0931 1428 char branchname[20];
fe4da5cc 1429
ddae0931 1430
1431 AliDetector::MakeBranch(option);
1432 sprintf(branchname,"%sCerenkov",GetName());
1433 if (fCerenkovs && gAlice->TreeH()) {
237c933d 1434 gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize);
ddae0931 1435 printf("Making Branch %s for Cerenkov Hits\n",branchname);
fe4da5cc 1436 }
ddae0931 1437
2e5f0f7b 1438 sprintf(branchname,"%sPadHits",GetName());
1439 if (fPadHits && gAlice->TreeH()) {
237c933d 1440 gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
2e5f0f7b 1441 printf("Making Branch %s for PadHits\n",branchname);
fe4da5cc 1442 }
1443
ddae0931 1444// one branch for digits per chamber
1445 Int_t i;
1446
237c933d 1447 for (i=0; i<kNCH ;i++) {
ddae0931 1448 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1449
1450 if (fDchambers && gAlice->TreeD()) {
237c933d 1451 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
ddae0931 1452 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
1453 }
fe4da5cc 1454 }
2e5f0f7b 1455
1456// one branch for raw clusters per chamber
237c933d 1457 for (i=0; i<kNCH ;i++) {
2e5f0f7b 1458 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1459
1460 if (fRawClusters && gAlice->TreeR()) {
237c933d 1461 gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
2e5f0f7b 1462 printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
1463 }
1464 }
1465
1466 // one branch for rec hits per chamber
237c933d 1467 for (i=0; i<kNCH ;i++) {
a4622d0f 1468 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1469
1470 if (fRecHits1D && gAlice->TreeR()) {
1471 gAlice->TreeR()->Branch(branchname,&((*fRecHits1D)[i]), kBufferSize);
1472 printf("Making Branch %s for 1D rec. hits in chamber %d\n",branchname,i+1);
1473 }
1474 }
1475 for (i=0; i<kNCH ;i++) {
1476 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
2e5f0f7b 1477
a4622d0f 1478 if (fRecHits3D && gAlice->TreeR()) {
1479 gAlice->TreeR()->Branch(branchname,&((*fRecHits3D)[i]), kBufferSize);
1480 printf("Making Branch %s for 3D rec. hits in chamber %d\n",branchname,i+1);
2e5f0f7b 1481 }
1482 }
a4622d0f 1483
fe4da5cc 1484}
1485
ddae0931 1486//___________________________________________
1487void AliRICH::SetTreeAddress()
fe4da5cc 1488{
237c933d 1489 // Set branch address for the Hits and Digits Tree.
2e5f0f7b 1490 char branchname[20];
1491 Int_t i;
1492
ddae0931 1493 AliDetector::SetTreeAddress();
1494
1495 TBranch *branch;
1496 TTree *treeH = gAlice->TreeH();
1497 TTree *treeD = gAlice->TreeD();
2e5f0f7b 1498 TTree *treeR = gAlice->TreeR();
ddae0931 1499
1500 if (treeH) {
2e5f0f7b 1501 if (fPadHits) {
1502 branch = treeH->GetBranch("RICHPadHits");
1503 if (branch) branch->SetAddress(&fPadHits);
fe4da5cc 1504 }
ddae0931 1505 if (fCerenkovs) {
1506 branch = treeH->GetBranch("RICHCerenkov");
1507 if (branch) branch->SetAddress(&fCerenkovs);
fe4da5cc 1508 }
ddae0931 1509 }
1510
1511 if (treeD) {
237c933d 1512 for (int i=0; i<kNCH; i++) {
ddae0931 1513 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1514 if (fDchambers) {
1515 branch = treeD->GetBranch(branchname);
1516 if (branch) branch->SetAddress(&((*fDchambers)[i]));
1517 }
fe4da5cc 1518 }
fe4da5cc 1519 }
2e5f0f7b 1520 if (treeR) {
237c933d 1521 for (i=0; i<kNCH; i++) {
2e5f0f7b 1522 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1523 if (fRawClusters) {
1524 branch = treeR->GetBranch(branchname);
1525 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1526 }
1527 }
1528
237c933d 1529 for (i=0; i<kNCH; i++) {
a4622d0f 1530 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1531 if (fRecHits1D) {
2e5f0f7b 1532 branch = treeR->GetBranch(branchname);
a4622d0f 1533 if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
2e5f0f7b 1534 }
1535 }
1536
a4622d0f 1537 for (i=0; i<kNCH; i++) {
1538 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1539 if (fRecHits3D) {
1540 branch = treeR->GetBranch(branchname);
1541 if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
1542 }
1543 }
1544
2e5f0f7b 1545 }
ddae0931 1546}
1547//___________________________________________
1548void AliRICH::ResetHits()
1549{
237c933d 1550 // Reset number of clusters and the cluster array for this detector
ddae0931 1551 AliDetector::ResetHits();
2e5f0f7b 1552 fNPadHits = 0;
1553 fNcerenkovs = 0;
1554 if (fPadHits) fPadHits->Clear();
ddae0931 1555 if (fCerenkovs) fCerenkovs->Clear();
fe4da5cc 1556}
1557
2e5f0f7b 1558
ddae0931 1559//____________________________________________
1560void AliRICH::ResetDigits()
1561{
237c933d 1562 //
1563 // Reset number of digits and the digits array for this detector
1564 //
1565 for ( int i=0;i<kNCH;i++ ) {
ddae0931 1566 if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
1567 if (fNdch) fNdch[i]=0;
1568 }
1569}
2e5f0f7b 1570
ddae0931 1571//____________________________________________
2e5f0f7b 1572void AliRICH::ResetRawClusters()
fe4da5cc 1573{
237c933d 1574 //
1575 // Reset number of raw clusters and the raw clust array for this detector
1576 //
1577 for ( int i=0;i<kNCH;i++ ) {
2e5f0f7b 1578 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
1579 if (fNrawch) fNrawch[i]=0;
fe4da5cc 1580 }
ddae0931 1581}
fe4da5cc 1582
2e5f0f7b 1583//____________________________________________
a4622d0f 1584void AliRICH::ResetRecHits1D()
fe4da5cc 1585{
237c933d 1586 //
1587 // Reset number of raw clusters and the raw clust array for this detector
1588 //
2e5f0f7b 1589
237c933d 1590 for ( int i=0;i<kNCH;i++ ) {
a4622d0f 1591 if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear();
1592 if (fNrechits1D) fNrechits1D[i]=0;
1593 }
1594}
1595
1596//____________________________________________
1597void AliRICH::ResetRecHits3D()
1598{
1599 //
1600 // Reset number of raw clusters and the raw clust array for this detector
1601 //
1602
1603 for ( int i=0;i<kNCH;i++ ) {
1604 if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear();
1605 if (fNrechits3D) fNrechits3D[i]=0;
2e5f0f7b 1606 }
fe4da5cc 1607}
1608
ddae0931 1609//___________________________________________
2e5f0f7b 1610void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
fe4da5cc 1611{
237c933d 1612
1613//
1614// Setter for the RICH geometry model
1615//
1616
1617
2e5f0f7b 1618 ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
fe4da5cc 1619}
1620
ddae0931 1621//___________________________________________
a2f7eaf6 1622void AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation)
ddae0931 1623{
237c933d 1624
1625//
1626// Setter for the RICH segmentation model
1627//
1628
a2f7eaf6 1629 ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation);
ddae0931 1630}
fe4da5cc 1631
ddae0931 1632//___________________________________________
2e5f0f7b 1633void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
fe4da5cc 1634{
237c933d 1635
1636//
1637// Setter for the RICH response model
1638//
1639
2e5f0f7b 1640 ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
fe4da5cc 1641}
1642
2e5f0f7b 1643void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
ddae0931 1644{
237c933d 1645
1646//
1647// Setter for the RICH reconstruction model (clusters)
1648//
1649
a2f7eaf6 1650 ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
fe4da5cc 1651}
1652
ddae0931 1653void AliRICH::SetNsec(Int_t id, Int_t nsec)
1654{
237c933d 1655
1656//
1657// Sets the number of padplanes
1658//
1659
2e5f0f7b 1660 ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
ddae0931 1661}
1662
1663
1664//___________________________________________
ddae0931 1665void AliRICH::StepManager()
fe4da5cc 1666{
237c933d 1667
452a64c6 1668// Full Step Manager
1669
1670 Int_t copy, id;
1671 static Int_t idvol;
1672 static Int_t vol[2];
1673 Int_t ipart;
1674 static Float_t hits[18];
1675 static Float_t ckovData[19];
1676 TLorentzVector position;
1677 TLorentzVector momentum;
1678 Float_t pos[3];
1679 Float_t mom[4];
1680 Float_t localPos[3];
1681 Float_t localMom[4];
1682 Float_t localTheta,localPhi;
1683 Float_t theta,phi;
1684 Float_t destep, step;
1685 Float_t ranf[2];
1686 Int_t nPads;
1687 Float_t coscerenkov;
1688 static Float_t eloss, xhit, yhit, tlength;
1689 const Float_t kBig=1.e10;
1690
1691 TClonesArray &lhits = *fHits;
452a64c6 1692 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
237c933d 1693
452a64c6 1694 //if (current->Energy()>1)
1695 //{
1696
1697 // Only gas gap inside chamber
1698 // Tag chambers and record hits when track enters
1699
1700 idvol=-1;
1701 id=gMC->CurrentVolID(copy);
1702 Float_t cherenkovLoss=0;
1703 //gAlice->KeepTrack(gAlice->CurrentTrack());
1704
1705 gMC->TrackPosition(position);
1706 pos[0]=position(0);
1707 pos[1]=position(1);
1708 pos[2]=position(2);
8cda28a5 1709 //bzero((char *)ckovData,sizeof(ckovData)*19);
452a64c6 1710 ckovData[1] = pos[0]; // X-position for hit
1711 ckovData[2] = pos[1]; // Y-position for hit
1712 ckovData[3] = pos[2]; // Z-position for hit
3cc8edee 1713 ckovData[6] = 0; // dummy track length
452a64c6 1714 //ckovData[11] = gAlice->CurrentTrack();
8cda28a5 1715
1716 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
452a64c6 1717
1718 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
1719
1720 /********************Store production parameters for Cerenkov photons************************/
1721//is it a Cerenkov photon?
8cda28a5 1722 if (gMC->TrackPid() == 50000050) {
452a64c6 1723
1724 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
1725 //{
1726 Float_t ckovEnergy = current->Energy();
1727 //energy interval for tracking
1728 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
1729 //if (ckovEnergy > 0)
1730 {
8cda28a5 1731 if (gMC->IsTrackEntering()){ //is track entering?
1732 //printf("Track entered (1)\n");
452a64c6 1733 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1734 { //is it in freo?
5cfcdf54 1735 if (gMC->IsNewTrack()){ //is it the first step?
8cda28a5 1736 //printf("I'm in!\n");
452a64c6 1737 Int_t mother = current->GetFirstMother();
1738
1739 //printf("Second Mother:%d\n",current->GetSecondMother());
1740
1741 ckovData[10] = mother;
1742 ckovData[11] = gAlice->CurrentTrack();
1743 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
8cda28a5 1744 //printf("Produced in FREO\n");
452a64c6 1745 fCkovNumber++;
1746 fFreonProd=1;
1747 //printf("Index: %d\n",fCkovNumber);
1748 } //first step question
1749 } //freo question
1750
5cfcdf54 1751 if (gMC->IsNewTrack()){ //is it first step?
452a64c6 1752 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
1753 {
1754 ckovData[12] = 2;
8cda28a5 1755 //printf("Produced in QUAR\n");
452a64c6 1756 } //quarz question
1757 } //first step question
1758
1759 //printf("Before %d\n",fFreonProd);
1760 } //track entering question
1761
1762 if (ckovData[12] == 1) //was it produced in Freon?
1763 //if (fFreonProd == 1)
1764 {
1765 if (gMC->IsTrackEntering()){ //is track entering?
8cda28a5 1766 //printf("Track entered (2)\n");
1767 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
1768 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
452a64c6 1769 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
1770 {
8cda28a5 1771 //printf("Got in META\n");
452a64c6 1772 gMC->TrackMomentum(momentum);
1773 mom[0]=momentum(0);
1774 mom[1]=momentum(1);
1775 mom[2]=momentum(2);
1776 mom[3]=momentum(3);
1777 // Z-position for hit
1778
1779
1780 /**************** Photons lost in second grid have to be calculated by hand************/
1781
1782 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1783 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
1784 gMC->Rndm(ranf, 1);
1785 //printf("grid calculation:%f\n",t);
1786 if (ranf[0] > t) {
5cfcdf54 1787 gMC->StopTrack();
452a64c6 1788 ckovData[13] = 5;
1789 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1790 //printf("Added One (1)!\n");
452a64c6 1791 //printf("Lost one in grid\n");
1792 }
1793 /**********************************************************************************/
1794 } //gap
1795
8cda28a5 1796 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
1797 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
452a64c6 1798 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
1799 {
8cda28a5 1800 //printf("Got in CSI\n");
452a64c6 1801 gMC->TrackMomentum(momentum);
1802 mom[0]=momentum(0);
1803 mom[1]=momentum(1);
1804 mom[2]=momentum(2);
1805 mom[3]=momentum(3);
1806
1807 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
1808 /***********************Cerenkov phtons (always polarised)*************************/
1809
1810 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1811 Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
1812 gMC->Rndm(ranf, 1);
1813 if (ranf[0] < t) {
5cfcdf54 1814 gMC->StopTrack();
452a64c6 1815 ckovData[13] = 6;
1816 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1817 //printf("Added One (2)!\n");
452a64c6 1818 //printf("Lost by Fresnel\n");
1819 }
1820 /**********************************************************************************/
1821 }
1822 } //track entering?
1823
1824
1825 /********************Evaluation of losses************************/
1826 /******************still in the old fashion**********************/
1827
5cfcdf54 1828 TArrayI procs;
1829 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
452a64c6 1830 for (Int_t i = 0; i < i1; ++i) {
1831 // Reflection loss
5cfcdf54 1832 if (procs[i] == kPLightReflection) { //was it reflected
452a64c6 1833 ckovData[13]=10;
1834 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1835 ckovData[13]=1;
1836 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1837 ckovData[13]=2;
5cfcdf54 1838 //gMC->StopTrack();
9dda3582 1839 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
452a64c6 1840 } //reflection question
8cda28a5 1841
452a64c6 1842 // Absorption loss
5cfcdf54 1843 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
8cda28a5 1844 //printf("Got in absorption\n");
452a64c6 1845 ckovData[13]=20;
1846 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1847 ckovData[13]=11;
1848 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1849 ckovData[13]=12;
1850 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
1851 ckovData[13]=13;
1852 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
1853 ckovData[13]=13;
1854
1855 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
1856 ckovData[13]=15;
1857
1858 // CsI inefficiency
1859 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
1860 ckovData[13]=16;
1861 }
5cfcdf54 1862 gMC->StopTrack();
452a64c6 1863 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1864 //printf("Added One (3)!\n");
452a64c6 1865 //printf("Added cerenkov %d\n",fCkovNumber);
1866 } //absorption question
1867
1868
1869 // Photon goes out of tracking scope
5cfcdf54 1870 else if (procs[i] == kPStop) { //is it below energy treshold?
452a64c6 1871 ckovData[13]=21;
5cfcdf54 1872 gMC->StopTrack();
452a64c6 1873 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1874 //printf("Added One (4)!\n");
452a64c6 1875 } // energy treshold question
1876 } //number of mechanisms cycle
1877 /**********************End of evaluation************************/
1878 } //freon production question
1879 } //energy interval question
1880 //}//inside the proximity gap question
1881 } //cerenkov photon question
1882
1883 /**************************************End of Production Parameters Storing*********************/
1884
1885
1886 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
1887
1888 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
1889 //printf("Cerenkov\n");
1890 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
1891 {
8cda28a5 1892 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
1893 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
1894 //printf("Got in CSI\n");
452a64c6 1895 if (gMC->Edep() > 0.){
1896 gMC->TrackPosition(position);
1897 gMC->TrackMomentum(momentum);
1898 pos[0]=position(0);
1899 pos[1]=position(1);
1900 pos[2]=position(2);
1901 mom[0]=momentum(0);
1902 mom[1]=momentum(1);
1903 mom[2]=momentum(2);
1904 mom[3]=momentum(3);
1905 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1906 Double_t rt = TMath::Sqrt(tc);
1907 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1908 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1909 gMC->Gmtod(pos,localPos,1);
1910 gMC->Gmtod(mom,localMom,2);
1911
1912 gMC->CurrentVolOffID(2,copy);
1913 vol[0]=copy;
1914 idvol=vol[0]-1;
1915
1916 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1917 //->Sector(localPos[0], localPos[2]);
1918 //printf("Sector:%d\n",sector);
1919
1920 /*if (gMC->TrackPid() == 50000051){
1921 fFeedbacks++;
1922 printf("Feedbacks:%d\n",fFeedbacks);
1923 }*/
1924
1925 ((AliRICHChamber*) (*fChambers)[idvol])
1926 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
1927 if(idvol<kNCH) {
1928 ckovData[0] = gMC->TrackPid(); // particle type
1929 ckovData[1] = pos[0]; // X-position for hit
1930 ckovData[2] = pos[1]; // Y-position for hit
1931 ckovData[3] = pos[2]; // Z-position for hit
1932 ckovData[4] = theta; // theta angle of incidence
1933 ckovData[5] = phi; // phi angle of incidence
1934 ckovData[8] = (Float_t) fNPadHits; // first padhit
1935 ckovData[9] = -1; // last pad hit
1936 ckovData[13] = 4; // photon was detected
1937 ckovData[14] = mom[0];
1938 ckovData[15] = mom[1];
1939 ckovData[16] = mom[2];
1940
1941 destep = gMC->Edep();
1942 gMC->SetMaxStep(kBig);
1943 cherenkovLoss += destep;
1944 ckovData[7]=cherenkovLoss;
1945
1946 nPads = MakePadHits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
1947 if (fNPadHits > (Int_t)ckovData[8]) {
1948 ckovData[8]= ckovData[8]+1;
1949 ckovData[9]= (Float_t) fNPadHits;
1950 }
1951
1952 ckovData[17] = nPads;
1953 //printf("nPads:%d",nPads);
1954
1955 //TClonesArray *Hits = RICH->Hits();
1956 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
1957 if (mipHit)
1958 {
1959 mom[0] = current->Px();
1960 mom[1] = current->Py();
1961 mom[2] = current->Pz();
1962 Float_t mipPx = mipHit->fMomX;
1963 Float_t mipPy = mipHit->fMomY;
1964 Float_t mipPz = mipHit->fMomZ;
1965
1966 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
1967 Float_t rt = TMath::Sqrt(r);
1968 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
1969 Float_t mipRt = TMath::Sqrt(mipR);
1970 if ((rt*mipRt) > 0)
1971 {
1972 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
1973 }
1974 else
1975 {
1976 coscerenkov = 0;
1977 }
1978 Float_t cherenkov = TMath::ACos(coscerenkov);
1979 ckovData[18]=cherenkov;
1980 }
1981 //if (sector != -1)
1982 //{
1983 AddHit(gAlice->CurrentTrack(),vol,ckovData);
1984 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1985 //printf("Added One (5)!\n");
452a64c6 1986 //}
1987 }
1988 }
1989 }
1990 }
1991
1992 /***********************************************End of photon hits*********************************************/
1993
1994
1995 /**********************************************Charged particles treatment*************************************/
1996
1997 else if (gMC->TrackCharge())
1998 //else if (1 == 1)
1999 {
2000//If MIP
2001 /*if (gMC->IsTrackEntering())
2002 {
2003 hits[13]=20;//is track entering?
2004 }*/
2005 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2006 {
2007 fFreonProd=1;
2008 }
2009
2010 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2011// Get current particle id (ipart), track position (pos) and momentum (mom)
2012
2013 gMC->CurrentVolOffID(3,copy);
2014 vol[0]=copy;
2015 idvol=vol[0]-1;
2016
2017 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2018 //->Sector(localPos[0], localPos[2]);
2019 //printf("Sector:%d\n",sector);
2020
2021 gMC->TrackPosition(position);
2022 gMC->TrackMomentum(momentum);
2023 pos[0]=position(0);
2024 pos[1]=position(1);
2025 pos[2]=position(2);
2026 mom[0]=momentum(0);
2027 mom[1]=momentum(1);
2028 mom[2]=momentum(2);
2029 mom[3]=momentum(3);
2030 gMC->Gmtod(pos,localPos,1);
2031 gMC->Gmtod(mom,localMom,2);
2032
2033 ipart = gMC->TrackPid();
2034 //
2035 // momentum loss and steplength in last step
2036 destep = gMC->Edep();
2037 step = gMC->TrackStep();
2038
2039 //
2040 // record hits when track enters ...
2041 if( gMC->IsTrackEntering()) {
2042// gMC->SetMaxStep(fMaxStepGas);
2043 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2044 Double_t rt = TMath::Sqrt(tc);
2045 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2046 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2047
2048
2049 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2050 Double_t localRt = TMath::Sqrt(localTc);
2051 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2052 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2053
2054 hits[0] = Float_t(ipart); // particle type
2055 hits[1] = localPos[0]; // X-position for hit
2056 hits[2] = localPos[1]; // Y-position for hit
2057 hits[3] = localPos[2]; // Z-position for hit
2058 hits[4] = localTheta; // theta angle of incidence
2059 hits[5] = localPhi; // phi angle of incidence
2060 hits[8] = (Float_t) fNPadHits; // first padhit
2061 hits[9] = -1; // last pad hit
2062 hits[13] = fFreonProd; // did id hit the freon?
2063 hits[14] = mom[0];
2064 hits[15] = mom[1];
2065 hits[16] = mom[2];
2066
2067 tlength = 0;
2068 eloss = 0;
2069 fFreonProd = 0;
2070
2071 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2072
2073
2074 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2075 xhit = localPos[0];
2076 yhit = localPos[2];
2077 // Only if not trigger chamber
2078 if(idvol<kNCH) {
2079 //
2080 // Initialize hit position (cursor) in the segmentation model
2081 ((AliRICHChamber*) (*fChambers)[idvol])
2082 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2083 }
2084 }
2085
2086 //
2087 // Calculate the charge induced on a pad (disintegration) in case
2088 //
2089 // Mip left chamber ...
2090 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2091 gMC->SetMaxStep(kBig);
2092 eloss += destep;
2093 tlength += step;
2094
2095
2096 // Only if not trigger chamber
2097 if(idvol<kNCH) {
2098 if (eloss > 0)
2099 {
2100 if(gMC->TrackPid() == kNeutron)
2101 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2102 nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip);
2103 hits[17] = nPads;
2104 //printf("nPads:%d",nPads);
2105 }
2106 }
2107
2108 hits[6]=tlength;
2109 hits[7]=eloss;
2110 if (fNPadHits > (Int_t)hits[8]) {
2111 hits[8]= hits[8]+1;
2112 hits[9]= (Float_t) fNPadHits;
2113 }
2114
2115 //if(sector !=-1)
2116 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2117 eloss = 0;
2118 //
2119 // Check additional signal generation conditions
2120 // defined by the segmentation
2121 // model (boundary crossing conditions)
2122 } else if
2123 (((AliRICHChamber*) (*fChambers)[idvol])
2124 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2125 {
2126 ((AliRICHChamber*) (*fChambers)[idvol])
2127 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2128 if (eloss > 0)
2129 {
2130 if(gMC->TrackPid() == kNeutron)
2131 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2132 nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip);
2133 hits[17] = nPads;
2134 //printf("Npads:%d",NPads);
2135 }
2136 xhit = localPos[0];
2137 yhit = localPos[2];
2138 eloss = destep;
2139 tlength += step ;
2140 //
2141 // nothing special happened, add up energy loss
2142 } else {
2143 eloss += destep;
2144 tlength += step ;
2145 }
2146 }
2147 }
2148 /*************************************************End of MIP treatment**************************************/
2149 //}
ddae0931 2150}
2e5f0f7b 2151
237c933d 2152void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
ddae0931 2153{
2e5f0f7b 2154
ddae0931 2155//
2156// Loop on chambers and on cathode planes
2157//
2e5f0f7b 2158 for (Int_t icat=1;icat<2;icat++) {
2159 gAlice->ResetDigits();
2160 gAlice->TreeD()->GetEvent(1); // spurious +1 ...
237c933d 2161 for (Int_t ich=0;ich<kNCH;ich++) {
2e5f0f7b 2162 AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
237c933d 2163 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
2164 if (pRICHdigits == 0)
2e5f0f7b 2165 continue;
2166 //
2167 // Get ready the current chamber stuff
2168 //
2169 AliRICHResponse* response = iChamber->GetResponseModel();
a2f7eaf6 2170 AliSegmentation* seg = iChamber->GetSegmentationModel();
2e5f0f7b 2171 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2172 if (seg) {
2173 rec->SetSegmentation(seg);
2174 rec->SetResponse(response);
237c933d 2175 rec->SetDigits(pRICHdigits);
2e5f0f7b 2176 rec->SetChamber(ich);
2177 if (nev==0) rec->CalibrateCOG();
2178 rec->FindRawClusters();
2179 }
2180 TClonesArray *fRch;
2181 fRch=RawClustAddress(ich);
2182 fRch->Sort();
2183 } // for ich
2184
2185 gAlice->TreeR()->Fill();
2186 TClonesArray *fRch;
237c933d 2187 for (int i=0;i<kNCH;i++) {
2e5f0f7b 2188 fRch=RawClustAddress(i);
2189 int nraw=fRch->GetEntriesFast();
2190 printf ("Chamber %d, raw clusters %d\n",i,nraw);
2191 }
2192
2193 ResetRawClusters();
2194
2195 } // for icat
2196
2197 char hname[30];
2198 sprintf(hname,"TreeR%d",nev);
33984590 2199 gAlice->TreeR()->Write(hname,kOverwrite,0);
2e5f0f7b 2200 gAlice->TreeR()->Reset();
2201
2202 //gObjectTable->Print();
ddae0931 2203}
2204
2205
2206//______________________________________________________________________________
2207void AliRICH::Streamer(TBuffer &R__b)
2208{
2209 // Stream an object of class AliRICH.
2e5f0f7b 2210 AliRICHChamber *iChamber;
a2f7eaf6 2211 AliSegmentation *segmentation;
2e5f0f7b 2212 AliRICHResponse *response;
ddae0931 2213 TClonesArray *digitsaddress;
2e5f0f7b 2214 TClonesArray *rawcladdress;
a4622d0f 2215 TClonesArray *rechitaddress1D;
2216 TClonesArray *rechitaddress3D;
ddae0931 2217
2218 if (R__b.IsReading()) {
2219 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2220 AliDetector::Streamer(R__b);
2e5f0f7b 2221 R__b >> fNPadHits;
2222 R__b >> fPadHits; // diff
2223 R__b >> fNcerenkovs;
2224 R__b >> fCerenkovs; // diff
ddae0931 2225 R__b >> fDchambers;
2e5f0f7b 2226 R__b >> fRawClusters;
a4622d0f 2227 R__b >> fRecHits1D; //diff
2228 R__b >> fRecHits3D; //diff
33984590 2229 R__b >> fDebugLevel; //diff
9dda3582 2230 R__b.ReadStaticArray(fNdch);
2231 R__b.ReadStaticArray(fNrawch);
a4622d0f 2232 R__b.ReadStaticArray(fNrechits1D);
2233 R__b.ReadStaticArray(fNrechits3D);
2e5f0f7b 2234//
ddae0931 2235 R__b >> fChambers;
2236// Stream chamber related information
237c933d 2237 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 2238 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 2239 iChamber->Streamer(R__b);
2e5f0f7b 2240 segmentation=iChamber->GetSegmentationModel();
2241 segmentation->Streamer(R__b);
2242 response=iChamber->GetResponseModel();
ddae0931 2243 response->Streamer(R__b);
2e5f0f7b 2244 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
2245 rawcladdress->Streamer(R__b);
a4622d0f 2246 rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i];
2247 rechitaddress1D->Streamer(R__b);
2248 rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i];
2249 rechitaddress3D->Streamer(R__b);
ddae0931 2250 digitsaddress=(TClonesArray*) (*fDchambers)[i];
2251 digitsaddress->Streamer(R__b);
2252 }
452a64c6 2253 R__b >> fDebugLevel;
2254 R__b >> fCkovNumber;
2255 R__b >> fCkovQuarz;
2256 R__b >> fCkovGap;
2257 R__b >> fCkovCsi;
2258 R__b >> fLostRfreo;
2259 R__b >> fLostRquar;
2260 R__b >> fLostAfreo;
2261 R__b >> fLostAquarz;
2262 R__b >> fLostAmeta;
2263 R__b >> fLostCsi;
2264 R__b >> fLostWires;
2265 R__b >> fFreonProd;
2266 R__b >> fMipx;
2267 R__b >> fMipy;
2268 R__b >> fFeedbacks;
2269 R__b >> fLostFresnel;
ddae0931 2270
2271 } else {
2272 R__b.WriteVersion(AliRICH::IsA());
2273 AliDetector::Streamer(R__b);
2e5f0f7b 2274 R__b << fNPadHits;
2275 R__b << fPadHits; // diff
2276 R__b << fNcerenkovs;
2277 R__b << fCerenkovs; // diff
ddae0931 2278 R__b << fDchambers;
2e5f0f7b 2279 R__b << fRawClusters;
a4622d0f 2280 R__b << fRecHits1D; //diff
2281 R__b << fRecHits3D; //diff
33984590 2282 R__b << fDebugLevel; //diff
237c933d 2283 R__b.WriteArray(fNdch, kNCH);
2284 R__b.WriteArray(fNrawch, kNCH);
a4622d0f 2285 R__b.WriteArray(fNrechits1D, kNCH);
2286 R__b.WriteArray(fNrechits3D, kNCH);
2287//
ddae0931 2288 R__b << fChambers;
2289// Stream chamber related information
237c933d 2290 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 2291 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 2292 iChamber->Streamer(R__b);
2e5f0f7b 2293 segmentation=iChamber->GetSegmentationModel();
2294 segmentation->Streamer(R__b);
2295 response=iChamber->GetResponseModel();
ddae0931 2296 response->Streamer(R__b);
2e5f0f7b 2297 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
2298 rawcladdress->Streamer(R__b);
a4622d0f 2299 rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i];
2300 rechitaddress1D->Streamer(R__b);
2301 rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i];
2302 rechitaddress3D->Streamer(R__b);
ddae0931 2303 digitsaddress=(TClonesArray*) (*fDchambers)[i];
2304 digitsaddress->Streamer(R__b);
2305 }
452a64c6 2306 R__b << fDebugLevel;
2307 R__b << fCkovNumber;
2308 R__b << fCkovQuarz;
2309 R__b << fCkovGap;
2310 R__b << fCkovCsi;
2311 R__b << fLostRfreo;
2312 R__b << fLostRquar;
2313 R__b << fLostAfreo;
2314 R__b << fLostAquarz;
2315 R__b << fLostAmeta;
2316 R__b << fLostCsi;
2317 R__b << fLostWires;
2318 R__b << fFreonProd;
2319 R__b << fMipx;
2320 R__b << fMipy;
2321 R__b << fFeedbacks;
2322 R__b << fLostFresnel;
fe4da5cc 2323 }
ddae0931 2324}
2e5f0f7b 2325AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
ddae0931 2326{
2327//
2328 // Initialise the pad iterator
2329 // Return the address of the first padhit for hit
2330 TClonesArray *theClusters = clusters;
2331 Int_t nclust = theClusters->GetEntriesFast();
2332 if (nclust && hit->fPHlast > 0) {
2e5f0f7b 2333 sMaxIterPad=Int_t(hit->fPHlast);
2334 sCurIterPad=Int_t(hit->fPHfirst);
2335 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2336 } else {
2337 return 0;
fe4da5cc 2338 }
ddae0931 2339
fe4da5cc 2340}
2341
2e5f0f7b 2342AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters)
fe4da5cc 2343{
237c933d 2344
2345 // Iterates over pads
2346
ddae0931 2347 sCurIterPad++;
2348 if (sCurIterPad <= sMaxIterPad) {
2e5f0f7b 2349 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2350 } else {
2351 return 0;
2352 }
fe4da5cc 2353}
2354
fe4da5cc 2355
d28b5fc2 2356void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
fe4da5cc 2357{
ddae0931 2358 // keep galice.root for signal and name differently the file for
2359 // background when add! otherwise the track info for signal will be lost !
2e5f0f7b 2360
c90dd3e2 2361 static Bool_t first=kTRUE;
237c933d 2362 static TFile *pFile;
2363 char *addBackground = strstr(option,"Add");
53a60579 2364 Int_t particle;
ddae0931 2365
2e5f0f7b 2366 FILE* points; //these will be the digits...
2367
2368 points=fopen("points.dat","w");
2369
2370 AliRICHChamber* iChamber;
a2f7eaf6 2371 AliSegmentation* segmentation;
ddae0931 2372
53a60579 2373 Int_t digitise=0;
ddae0931 2374 Int_t trk[50];
2375 Int_t chtrk[50];
2376 TObjArray *list=new TObjArray;
237c933d 2377 static TClonesArray *pAddress=0;
2378 if(!pAddress) pAddress=new TClonesArray("TVector",1000);
2e5f0f7b 2379 Int_t digits[5];
ddae0931 2380
237c933d 2381 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
a2f7eaf6 2382 AliHitMap* pHitMap[10];
2e5f0f7b 2383 Int_t i;
237c933d 2384 for (i=0; i<10; i++) {pHitMap[i]=0;}
2385 if (addBackground ) {
ddae0931 2386 if(first) {
2387 fFileName=filename;
2388 cout<<"filename"<<fFileName<<endl;
237c933d 2389 pFile=new TFile(fFileName);
ddae0931 2390 cout<<"I have opened "<<fFileName<<" file "<<endl;
2e5f0f7b 2391 fHits2 = new TClonesArray("AliRICHHit",1000 );
2392 fClusters2 = new TClonesArray("AliRICHPadHit",10000);
c90dd3e2 2393 first=kFALSE;
ddae0931 2394 }
237c933d 2395 pFile->cd();
ddae0931 2396 // Get Hits Tree header from file
2397 if(fHits2) fHits2->Clear();
2398 if(fClusters2) fClusters2->Clear();
2e5f0f7b 2399 if(TrH1) delete TrH1;
2400 TrH1=0;
2401
ddae0931 2402 char treeName[20];
2403 sprintf(treeName,"TreeH%d",nev);
2e5f0f7b 2404 TrH1 = (TTree*)gDirectory->Get(treeName);
2405 if (!TrH1) {
ddae0931 2406 printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
2407 }
2408 // Set branch addresses
2409 TBranch *branch;
2410 char branchname[20];
2411 sprintf(branchname,"%s",GetName());
2e5f0f7b 2412 if (TrH1 && fHits2) {
2413 branch = TrH1->GetBranch(branchname);
ddae0931 2414 if (branch) branch->SetAddress(&fHits2);
2415 }
2e5f0f7b 2416 if (TrH1 && fClusters2) {
2417 branch = TrH1->GetBranch("RICHCluster");
ddae0931 2418 if (branch) branch->SetAddress(&fClusters2);
2419 }
2420 }
33984590 2421
a2f7eaf6 2422 AliHitMap* hm;
2e5f0f7b 2423 Int_t countadr=0;
33984590 2424 Int_t counter=0;
2425 for (i =0; i<kNCH; i++) {
2426 iChamber=(AliRICHChamber*) (*fChambers)[i];
2427 segmentation=iChamber->GetSegmentationModel(1);
2428 pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
2429 }
2430 //
2431 // Loop over tracks
2432 //
2433
2434 TTree *treeH = gAlice->TreeH();
2435 Int_t ntracks =(Int_t) treeH->GetEntries();
2436 for (Int_t track=0; track<ntracks; track++) {
2437 gAlice->ResetHits();
2438 treeH->GetEvent(track);
2439 //
2440 // Loop over hits
2441 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2442 mHit;
2443 mHit=(AliRICHHit*)pRICH->NextHit())
2444 {
2445
33984590 2446 Int_t nch = mHit->fChamber-1; // chamber number
53a60579 2447 Int_t index = mHit->Track();
33984590 2448 if (nch >kNCH) continue;
2449 iChamber = &(pRICH->Chamber(nch));
2450
53a60579 2451 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
33984590 2452
53a60579 2453 if (current->GetPdgCode() >= 50000050)
2454 {
2455 TParticle *motherofcurrent = (TParticle*)(*gAlice->Particles())[current->GetFirstMother()];
2456 particle = motherofcurrent->GetPdgCode();
2457 }
2458 else
2459 {
2460 particle = current->GetPdgCode();
2461 }
2462
33984590 2463
2464 //printf("Flag:%d\n",flag);
2465 //printf("Track:%d\n",track);
2466 //printf("Particle:%d\n",particle);
2467
53a60579 2468 digitise=1;
33984590 2469
2470 if (flag == 1)
53a60579 2471 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2472 digitise=0;
33984590 2473
2474 if (flag == 2)
2475 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2476 || TMath::Abs(particle)==311)
53a60579 2477 digitise=0;
33984590 2478
2479 if (flag == 3 && TMath::Abs(particle)==2212)
53a60579 2480 digitise=0;
33984590 2481
2482 if (flag == 4 && TMath::Abs(particle)==13)
53a60579 2483 digitise=0;
33984590 2484
2485 if (flag == 5 && TMath::Abs(particle)==11)
53a60579 2486 digitise=0;
33984590 2487
2488 if (flag == 6 && TMath::Abs(particle)==2112)
53a60579 2489 digitise=0;
33984590 2490
2491
53a60579 2492 //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise);
33984590 2493
2494
53a60579 2495 if (digitise)
ddae0931 2496 {
d28b5fc2 2497
33984590 2498 //
2499 // Loop over pad hits
2500 for (AliRICHPadHit* mPad=
2501 (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
2502 mPad;
2503 mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
d28b5fc2 2504 {
33984590 2505 Int_t ipx = mPad->fPadX; // pad number on X
2506 Int_t ipy = mPad->fPadY; // pad number on Y
2507 Int_t iqpad = mPad->fQpad; // charge per pad
2508 //
2509 //
2510 //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
2511
a2f7eaf6 2512 Float_t thex, they, thez;
176a9917 2513 segmentation=iChamber->GetSegmentationModel(0);
a2f7eaf6 2514 segmentation->GetPadC(ipx,ipy,thex,they,thez);
33984590 2515 new((*pAddress)[countadr++]) TVector(2);
2516 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
2517 trinfo(0)=(Float_t)track;
2518 trinfo(1)=(Float_t)iqpad;
2519
2520 digits[0]=ipx;
2521 digits[1]=ipy;
2522 digits[2]=iqpad;
2523
2524 AliRICHTransientDigit* pdigit;
2525 // build the list of fired pads and update the info
2526 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
2527 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
2528 pHitMap[nch]->SetHit(ipx, ipy, counter);
2529 counter++;
2530 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
2531 // list of tracks
2532 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
2533 trlist->Add(&trinfo);
2534 } else {
2535 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
2536 // update charge
2537 (*pdigit).fSignal+=iqpad;
2538 // update list of tracks
2539 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
2540 Int_t lastEntry=trlist->GetLast();
2541 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
2542 TVector &ptrk=*ptrkP;
2543 Int_t lastTrack=Int_t(ptrk(0));
2544 Int_t lastCharge=Int_t(ptrk(1));
2545 if (lastTrack==track) {
2546 lastCharge+=iqpad;
2547 trlist->RemoveAt(lastEntry);
2548 trinfo(0)=lastTrack;
2549 trinfo(1)=lastCharge;
2550 trlist->AddAt(&trinfo,lastEntry);
2551 } else {
2552 trlist->Add(&trinfo);
2553 }
2554 // check the track list
2555 Int_t nptracks=trlist->GetEntriesFast();
2556 if (nptracks > 2) {
2557 printf("Attention - tracks: %d (>2)\n",nptracks);
2558 //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
2559 for (Int_t tr=0;tr<nptracks;tr++) {
2560 TVector *pptrkP=(TVector*)trlist->At(tr);
2561 TVector &pptrk=*pptrkP;
2562 trk[tr]=Int_t(pptrk(0));
2563 chtrk[tr]=Int_t(pptrk(1));
2564 }
2565 } // end if nptracks
2566 } // end if pdigit
2567 } //end loop over clusters
2568 }// track type condition
2569 } // hit loop
2570 } // track loop
2571
2572 // open the file with background
2573
2574 if (addBackground ) {
2575 ntracks =(Int_t)TrH1->GetEntries();
2576 //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
2577 //printf("background - Start loop over tracks \n");
2578 //
2579 // Loop over tracks
2580 //
2581 for (Int_t trak=0; trak<ntracks; trak++) {
2582 if (fHits2) fHits2->Clear();
2583 if (fClusters2) fClusters2->Clear();
2584 TrH1->GetEvent(trak);
2585 //
2586 // Loop over hits
2587 AliRICHHit* mHit;
2588 for(int j=0;j<fHits2->GetEntriesFast();++j)
2589 {
2590 mHit=(AliRICHHit*) (*fHits2)[j];
2591 Int_t nch = mHit->fChamber-1; // chamber number
2592 if (nch >6) continue;
2593 iChamber = &(pRICH->Chamber(nch));
2594 Int_t rmin = (Int_t)iChamber->RInner();
2595 Int_t rmax = (Int_t)iChamber->ROuter();
2596 //
2597 // Loop over pad hits
2598 for (AliRICHPadHit* mPad=
2599 (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
2600 mPad;
2601 mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2))
2602 {
33984590 2603 Int_t ipx = mPad->fPadX; // pad number on X
2604 Int_t ipy = mPad->fPadY; // pad number on Y
2605 Int_t iqpad = mPad->fQpad; // charge per pad
2e5f0f7b 2606
a2f7eaf6 2607 Float_t thex, they, thez;
176a9917 2608 segmentation=iChamber->GetSegmentationModel(0);
a2f7eaf6 2609 segmentation->GetPadC(ipx,ipy,thex,they,thez);
33984590 2610 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
2611 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
2612 new((*pAddress)[countadr++]) TVector(2);
2613 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
2614 trinfo(0)=-1; // tag background
2615 trinfo(1)=-1;
2616 digits[0]=ipx;
2617 digits[1]=ipy;
2618 digits[2]=iqpad;
2619 if (trak <4 && nch==0)
2620 printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
2621 pHitMap[nch]->TestHit(ipx, ipy),trak);
2622 AliRICHTransientDigit* pdigit;
2623 // build the list of fired pads and update the info
2624 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
2625 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
2626
2627 pHitMap[nch]->SetHit(ipx, ipy, counter);
2628 counter++;
2629 printf("bgr new elem in list - counter %d\n",counter);
2630
2631 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
2632 // list of tracks
2633 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
2634 trlist->Add(&trinfo);
2635 } else {
2636 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
2637 // update charge
2638 (*pdigit).fSignal+=iqpad;
2639 // update list of tracks
2640 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
2641 Int_t lastEntry=trlist->GetLast();
2642 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
2643 TVector &ptrk=*ptrkP;
2644 Int_t lastTrack=Int_t(ptrk(0));
2645 if (lastTrack==-1) {
2646 continue;
2647 } else {
2648 trlist->Add(&trinfo);
2649 }
2650 // check the track list
2651 Int_t nptracks=trlist->GetEntriesFast();
2652 if (nptracks > 0) {
2653 for (Int_t tr=0;tr<nptracks;tr++) {
2654 TVector *pptrkP=(TVector*)trlist->At(tr);
2655 TVector &pptrk=*pptrkP;
2656 trk[tr]=Int_t(pptrk(0));
2657 chtrk[tr]=Int_t(pptrk(1));
2658 }
2659 } // end if nptracks
2660 } // end if pdigit
2661 } //end loop over clusters
2662 } // hit loop
2663 } // track loop
ddae0931 2664 TTree *fAli=gAlice->TreeK();
237c933d 2665 if (fAli) pFile =fAli->GetCurrentFile();
2666 pFile->cd();
33984590 2667 } // if Add
2668
2669 Int_t tracks[10];
2670 Int_t charges[10];
2671 //cout<<"Start filling digits \n "<<endl;
2672 Int_t nentries=list->GetEntriesFast();
2673 //printf(" \n \n nentries %d \n",nentries);
2674
2675 // start filling the digits
2676
2677 for (Int_t nent=0;nent<nentries;nent++) {
2678 AliRICHTransientDigit *address=(AliRICHTransientDigit*)list->At(nent);
2679 if (address==0) continue;
2680
2681 Int_t ich=address->fChamber;
2682 Int_t q=address->fSignal;
2683 iChamber=(AliRICHChamber*) (*fChambers)[ich];
2684 AliRICHResponse * response=iChamber->GetResponseModel();
2685 Int_t adcmax= (Int_t) response->MaxAdc();
2686
2687
2688 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
2689 //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY));
2690 Int_t pedestal = iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY);
ddae0931 2691
33984590 2692 //printf("Pedestal:%d\n",pedestal);
2693 //Int_t pedestal=0;
9dda3582 2694 Float_t treshold = (pedestal + 4*2.2);
33984590 2695
9dda3582 2696 Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
33984590 2697 Float_t noise = gRandom->Gaus(0, meanNoise);
2698 q+=(Int_t)(noise + pedestal);
2699 //q+=(Int_t)(noise);
2700 // magic number to be parametrised !!!
2701 if ( q <= treshold)
2702 {
2703 q = q - pedestal;
2704 continue;
2e5f0f7b 2705 }
33984590 2706 q = q - pedestal;
2707 if ( q >= adcmax) q=adcmax;
2708 digits[0]=address->fPadX;
2709 digits[1]=address->fPadY;
2710 digits[2]=q;
2711
2712 TObjArray* trlist=(TObjArray*)address->TrackList();
2713 Int_t nptracks=trlist->GetEntriesFast();
2714
2715 // this was changed to accomodate the real number of tracks
2716 if (nptracks > 10) {
2717 cout<<"Attention - tracks > 10 "<<nptracks<<endl;
2718 nptracks=10;
2719 }
2720 if (nptracks > 2) {
2721 printf("Attention - tracks > 2 %d \n",nptracks);
2722 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
2723 //icat,ich,digits[0],digits[1],q);
2724 }
2725 for (Int_t tr=0;tr<nptracks;tr++) {
2726 TVector *ppP=(TVector*)trlist->At(tr);
2727 TVector &pp =*ppP;
2728 tracks[tr]=Int_t(pp(0));
2729 charges[tr]=Int_t(pp(1));
2730 } //end loop over list of tracks for one pad
2731 if (nptracks < 10 ) {
2732 for (Int_t t=nptracks; t<10; t++) {
2733 tracks[t]=0;
2734 charges[t]=0;
ddae0931 2735 }
33984590 2736 }
2737 //write file
2738 if (ich==2)
2739 fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
2740
2741 // fill digits
2742 pRICH->AddDigits(ich,tracks,charges,digits);
2743 }
2744 gAlice->TreeD()->Fill();
2745
2746 list->Delete();
2747 for(Int_t ii=0;ii<kNCH;++ii) {
2748 if (pHitMap[ii]) {
2749 hm=pHitMap[ii];
2750 delete hm;
2751 pHitMap[ii]=0;
2752 }
2753 }
2754
2755 //TTree *TD=gAlice->TreeD();
2756 //Stat_t ndig=TD->GetEntries();
2757 //cout<<"number of digits "<<ndig<<endl;
2758 TClonesArray *fDch;
2759 for (int k=0;k<kNCH;k++) {
2760 fDch= pRICH->DigitsAddress(k);
2761 int ndigit=fDch->GetEntriesFast();
2762 printf ("Chamber %d digits %d \n",k,ndigit);
2763 }
2764 pRICH->ResetDigits();
ddae0931 2765 char hname[30];
2766 sprintf(hname,"TreeD%d",nev);
33984590 2767 gAlice->TreeD()->Write(hname,kOverwrite,0);
2768
2769 // reset tree
2770 // gAlice->TreeD()->Reset();
2e5f0f7b 2771 delete list;
237c933d 2772 pAddress->Clear();
33984590 2773 // gObjectTable->Print();
2e5f0f7b 2774}
ddae0931 2775
237c933d 2776AliRICH& AliRICH::operator=(const AliRICH& rhs)
2777{
2778// Assignment operator
2779 return *this;
2780
2781}
ddae0931 2782
237c933d 2783Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
2784{
2785//
2786// Calls the charge disintegration method of the current chamber and adds
2787// the simulated cluster to the root treee
2788//
176a9917 2789 Int_t clhits[5];
2790 Float_t newclust[4][500];
237c933d 2791 Int_t nnew;
2792
2793//
2794// Integrated pulse height on chamber
2795
2796 clhits[0]=fNhits+1;
2797
2798 ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
2799 Int_t ic=0;
2800
2801//
2802// Add new clusters
2803 for (Int_t i=0; i<nnew; i++) {
176a9917 2804 if (Int_t(newclust[0][i]) > 0) {
237c933d 2805 ic++;
237c933d 2806// Cluster Charge
176a9917 2807 clhits[1] = Int_t(newclust[0][i]);
237c933d 2808// Pad: ix
176a9917 2809 clhits[2] = Int_t(newclust[1][i]);
237c933d 2810// Pad: iy
176a9917 2811 clhits[3] = Int_t(newclust[2][i]);
237c933d 2812// Pad: chamber sector
176a9917 2813 clhits[4] = Int_t(newclust[3][i]);
237c933d 2814
2815 AddPadHit(clhits);
2816 }
2817 }
2818return nnew;
2819}
ddae0931 2820