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