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