1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.53 2001/08/30 09:51:23 hristov
19 The operator[] is replaced by At() or AddAt() in case of TObjArray.
21 Revision 1.52 2001/05/16 14:57:20 alibrary
22 New files for folders and Stack
24 Revision 1.51 2001/05/14 10:18:55 hristov
25 Default arguments declared once
27 Revision 1.50 2001/05/10 14:44:16 jbarbosa
28 Corrected some overlaps (thanks I. Hrivnacovna).
30 Revision 1.49 2001/05/10 12:23:49 jbarbosa
31 Repositioned the RICH modules.
32 Eliminated magic numbers.
33 Incorporated diagnostics (from macros).
35 Revision 1.48 2001/03/15 10:35:00 jbarbosa
36 Corrected bug in MakeBranch (was using a different version of STEER)
38 Revision 1.47 2001/03/14 18:13:56 jbarbosa
39 Several changes to adapt to new IO.
40 Removed digitising function, using AliRICHMerger::Digitise from now on.
42 Revision 1.46 2001/03/12 17:46:33 hristov
43 Changes needed on Sun with CC 5.0
45 Revision 1.45 2001/02/27 22:11:46 jbarbosa
46 Testing TreeS, removing of output.
48 Revision 1.44 2001/02/27 15:19:12 jbarbosa
49 Transition to SDigits.
51 Revision 1.43 2001/02/23 17:19:06 jbarbosa
52 Corrected photocathode definition in BuildGeometry().
54 Revision 1.42 2001/02/13 20:07:23 jbarbosa
55 Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum
56 when entering the freon. Corrected calls to particle stack.
58 Revision 1.41 2001/01/26 20:00:20 hristov
59 Major upgrade of AliRoot code
61 Revision 1.40 2001/01/24 20:58:03 jbarbosa
62 Enhanced BuildGeometry. Now the photocathodes are drawn.
64 Revision 1.39 2001/01/22 21:40:24 jbarbosa
65 Removing magic numbers
67 Revision 1.37 2000/12/20 14:07:25 jbarbosa
68 Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
70 Revision 1.36 2000/12/18 17:45:54 jbarbosa
71 Cleaned up PadHits object.
73 Revision 1.35 2000/12/15 16:49:40 jbarbosa
74 Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
76 Revision 1.34 2000/11/10 18:12:12 jbarbosa
77 Bug fix for AliRICHCerenkov (thanks to P. Hristov)
79 Revision 1.33 2000/11/02 10:09:01 jbarbosa
80 Minor bug correction (some pointers were not initialised in the default constructor)
82 Revision 1.32 2000/11/01 15:32:55 jbarbosa
83 Updated to handle both reconstruction algorithms.
85 Revision 1.31 2000/10/26 20:18:33 jbarbosa
86 Supports for methane and freon vessels
88 Revision 1.30 2000/10/24 13:19:12 jbarbosa
91 Revision 1.29 2000/10/19 19:39:25 jbarbosa
92 Some more changes to geometry. Further correction of digitisation "per part. type"
94 Revision 1.28 2000/10/17 20:50:57 jbarbosa
95 Inversed digtise by particle type (now, only the selected particle type is not digitsed).
96 Corrected several geometry minor bugs.
97 Added new parameter (opaque quartz thickness).
99 Revision 1.27 2000/10/11 10:33:55 jbarbosa
100 Corrected bug introduced by earlier revisions (CerenkovData array cannot be reset to zero on wach call of StepManager)
102 Revision 1.26 2000/10/03 21:44:08 morsch
103 Use AliSegmentation and AliHit abstract base classes.
105 Revision 1.25 2000/10/02 21:28:12 fca
106 Removal of useless dependecies via forward declarations
108 Revision 1.24 2000/10/02 15:43:17 jbarbosa
109 Fixed forward declarations.
110 Fixed honeycomb density.
111 Fixed cerenkov storing.
114 Revision 1.23 2000/09/13 10:42:14 hristov
115 Minor corrections for HP, DEC and Sun; strings.h included
117 Revision 1.22 2000/09/12 18:11:13 fca
118 zero hits area before using
120 Revision 1.21 2000/07/21 10:21:07 morsch
121 fNrawch = 0; and fNrechits = 0; in the default constructor.
123 Revision 1.20 2000/07/10 15:28:39 fca
124 Correction of the inheritance scheme
126 Revision 1.19 2000/06/30 16:29:51 dibari
127 Added kDebugLevel variable to control output size on demand
129 Revision 1.18 2000/06/12 15:15:46 jbarbosa
132 Revision 1.17 2000/06/09 14:58:37 jbarbosa
133 New digitisation per particle type
135 Revision 1.16 2000/04/19 12:55:43 morsch
136 Newly structured and updated version (JB, AM)
141 ////////////////////////////////////////////////
142 // Manager and hits classes for set:RICH //
143 ////////////////////////////////////////////////
151 #include <TObjArray.h>
154 #include <TParticle.h>
155 #include <TGeometry.h>
163 #include <iostream.h>
167 #include "AliSegmentation.h"
168 #include "AliRICHSegmentationV0.h"
169 #include "AliRICHHit.h"
170 #include "AliRICHCerenkov.h"
171 #include "AliRICHSDigit.h"
172 #include "AliRICHDigit.h"
173 #include "AliRICHTransientDigit.h"
174 #include "AliRICHRawCluster.h"
175 #include "AliRICHRecHit1D.h"
176 #include "AliRICHRecHit3D.h"
177 #include "AliRICHHitMapA1.h"
178 #include "AliRICHClusterFinder.h"
179 #include "AliRICHMerger.h"
183 #include "AliConst.h"
185 #include "AliPoints.h"
186 #include "AliCallf77.h"
189 // Static variables for the pad-hit iterator routines
190 static Int_t sMaxIterPad=0;
191 static Int_t sCurIterPad=0;
195 //___________________________________________
198 // Default constructor for RICH manager class
211 for (Int_t i=0; i<7; i++)
223 //___________________________________________
224 AliRICH::AliRICH(const char *name, const char *title)
225 : AliDetector(name,title)
229 <img src="gif/alirich.gif">
233 fHits = new TClonesArray("AliRICHHit",1000 );
234 gAlice->AddHitList(fHits);
235 fSDigits = new TClonesArray("AliRICHSDigit",100000);
236 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
237 gAlice->AddHitList(fCerenkovs);
238 //gAlice->AddHitList(fHits);
243 //fNdch = new Int_t[kNCH];
245 fDchambers = new TObjArray(kNCH);
247 fRecHits1D = new TObjArray(kNCH);
248 fRecHits3D = new TObjArray(kNCH);
252 for (i=0; i<kNCH ;i++) {
253 //PH (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
254 fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i);
258 //fNrawch = new Int_t[kNCH];
260 fRawClusters = new TObjArray(kNCH);
261 //printf("Created fRwClusters with adress:%p",fRawClusters);
263 for (i=0; i<kNCH ;i++) {
264 //PH (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
265 fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i);
269 //fNrechits = new Int_t[kNCH];
271 for (i=0; i<kNCH ;i++) {
272 //PH (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
273 fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
275 for (i=0; i<kNCH ;i++) {
276 //PH (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
277 fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
279 //printf("Created fRecHits with adress:%p",fRecHits);
282 SetMarkerColor(kRed);
284 /*fChambers = new TObjArray(kNCH);
285 for (i=0; i<kNCH; i++)
286 (*fChambers)[i] = new AliRICHChamber();*/
291 AliRICH::AliRICH(const AliRICH& RICH)
297 //___________________________________________
301 // Destructor of RICH manager class
308 //PH Delete TObjArrays
314 fDchambers->Delete();
318 fRawClusters->Delete();
322 fRecHits1D->Delete();
326 fRecHits3D->Delete();
333 //_____________________________________________________________________________
334 Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
337 // Calls the charge disintegration method of the current chamber and adds
338 // the simulated cluster to the root treee
341 Float_t newclust[4][500];
345 // Integrated pulse height on chamber
349 //PH ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
350 ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
355 for (Int_t i=0; i<nnew; i++) {
356 if (Int_t(newclust[0][i]) > 0) {
359 clhits[1] = Int_t(newclust[0][i]);
361 clhits[2] = Int_t(newclust[1][i]);
363 clhits[3] = Int_t(newclust[2][i]);
364 // Pad: chamber sector
365 clhits[4] = Int_t(newclust[3][i]);
367 //printf(" %d %d %d %d %d\n", clhits[0], clhits[1], clhits[2], clhits[3], clhits[4]);
375 gAlice->TreeS()->Fill();
376 gAlice->TreeS()->Write(0,TObject::kOverwrite);
377 //printf("Filled SDigits...\n");
382 //___________________________________________
383 void AliRICH::Hits2SDigits()
386 // Dummy: sdigits are created during transport.
387 // Called from alirun.
389 int nparticles = gAlice->GetNtrack();
390 cout << "Particles (RICH):" <<nparticles<<endl;
391 if (nparticles > 0) printf("SDigits were already generated.\n");
395 //___________________________________________
396 void AliRICH::SDigits2Digits(Int_t nev, Int_t flag)
401 // Called from macro. Multiple events, more functionality.
403 AliRICHChamber* iChamber;
405 printf("Generating tresholds...\n");
407 for(Int_t i=0;i<7;i++) {
408 iChamber = &(Chamber(i));
409 iChamber->GenerateTresholds();
412 int nparticles = gAlice->GetNtrack();
417 fMerger->Digitise(nev,flag);
420 //Digitise(nev,flag);
422 //___________________________________________
423 void AliRICH::SDigits2Digits()
428 // Called from alirun, single event only.
430 AliRICHChamber* iChamber;
432 printf("Generating tresholds...\n");
434 for(Int_t i=0;i<7;i++) {
435 iChamber = &(Chamber(i));
436 iChamber->GenerateTresholds();
439 int nparticles = gAlice->GetNtrack();
440 cout << "Particles (RICH):" <<nparticles<<endl;
445 fMerger->Digitise(0,0);
449 //___________________________________________
450 void AliRICH::Digits2Reco()
454 // Called from alirun, single event only.
456 int nparticles = gAlice->GetNtrack();
457 cout << "Particles (RICH):" <<nparticles<<endl;
458 if (nparticles > 0) FindClusters(0,0);
462 //___________________________________________
463 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
467 // Adds a hit to the Hits list
469 TClonesArray &lhits = *fHits;
470 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
472 //_____________________________________________________________________________
473 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
477 // Adds a RICH cerenkov hit to the Cerenkov Hits list
480 TClonesArray &lcerenkovs = *fCerenkovs;
481 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
482 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
484 //___________________________________________
485 void AliRICH::AddSDigit(Int_t *clhits)
489 // Add a RICH pad hit to the list
492 //printf("fsdigits:%p, data: %d\n",fSDigits,clhits[2]);
493 TClonesArray &lSDigits = *fSDigits;
494 new(lSDigits[fNSDigits++]) AliRICHSDigit(clhits);
496 //_____________________________________________________________________________
497 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
501 // Add a RICH digit to the list
504 //printf("fdigits:%p, data: %d\n",((TClonesArray*)(*fDchambers)[id]),digits[0]);
505 //PH TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
506 TClonesArray &ldigits = *((TClonesArray*)fDchambers->At(id));
507 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
510 //_____________________________________________________________________________
511 void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
514 // Add a RICH digit to the list
517 //PH TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
518 TClonesArray &lrawcl = *((TClonesArray*)fRawClusters->At(id));
519 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
522 //_____________________________________________________________________________
523 void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
527 // Add a RICH reconstructed hit to the list
530 //PH TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
531 TClonesArray &lrec1D = *((TClonesArray*)fRecHits1D->At(id));
532 new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
535 //_____________________________________________________________________________
536 void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
540 // Add a RICH reconstructed hit to the list
543 //PH TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
544 TClonesArray &lrec3D = *((TClonesArray*)fRecHits3D->At(id));
545 new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
548 //___________________________________________
549 void AliRICH::BuildGeometry()
554 // Builds a TNode geometry for event display
556 TNode *node, *subnode, *top;
558 const int kColorRICH = kRed;
560 top=gAlice->GetGeometry()->GetNode("alice");
562 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
563 AliRICHSegmentationV0* segmentation;
564 AliRICHChamber* iChamber;
565 AliRICHGeometry* geometry;
567 iChamber = &(pRICH->Chamber(0));
568 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
569 geometry=iChamber->GetGeometryModel();
571 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
573 Float_t padplane_width = segmentation->GetPadPlaneWidth();
574 Float_t padplane_length = segmentation->GetPadPlaneLength();
576 //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());
578 new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2);
580 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2);
581 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength());
583 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
584 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
585 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
586 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
587 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
588 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
589 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
591 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
593 new TRotMatrix("rot993","rot993",90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
594 new TRotMatrix("rot994","rot994",90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
595 new TRotMatrix("rot995","rot995",90., 0. , 90. , 90. , 0. , 0. );
596 new TRotMatrix("rot996","rot996",90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
597 new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
598 new TRotMatrix("rot998","rot998",90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
599 new TRotMatrix("rot999","rot999",90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
601 Float_t pos1[3]={0. , offset*cosphi , offset*sinphi};
602 Float_t pos2[3]={offset*sintheta , offset*costheta , 0. };
603 Float_t pos3[3]={0. , offset , 0.};
604 Float_t pos4[3]={-offset*sintheta , offset*costheta , 0.};
605 Float_t pos5[3]={offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
606 Float_t pos6[3]={0. , offset*cosphi , -offset*sinphi};
607 Float_t pos7[3]={ -offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
611 //Float_t pos1[3]={0,471.8999,165.2599};
612 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
613 //new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
614 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
615 node->SetLineColor(kColorRICH);
617 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
618 subnode->SetLineColor(kGreen);
619 fNodes->Add(subnode);
620 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,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(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
624 subnode->SetLineColor(kGreen);
625 fNodes->Add(subnode);
626 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
627 subnode->SetLineColor(kGreen);
628 fNodes->Add(subnode);
629 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
630 subnode->SetLineColor(kGreen);
631 fNodes->Add(subnode);
632 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
633 subnode->SetLineColor(kGreen);
634 fNodes->Add(subnode);
639 //Float_t pos2[3]={171,470,0};
640 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
641 //new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
642 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
643 node->SetLineColor(kColorRICH);
645 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
646 subnode->SetLineColor(kGreen);
647 fNodes->Add(subnode);
648 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
649 subnode->SetLineColor(kGreen);
650 fNodes->Add(subnode);
651 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
652 subnode->SetLineColor(kGreen);
653 fNodes->Add(subnode);
654 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
655 subnode->SetLineColor(kGreen);
656 fNodes->Add(subnode);
657 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
658 subnode->SetLineColor(kGreen);
659 fNodes->Add(subnode);
660 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
661 subnode->SetLineColor(kGreen);
662 fNodes->Add(subnode);
667 //Float_t pos3[3]={0,500,0};
668 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
669 //new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
670 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
671 node->SetLineColor(kColorRICH);
673 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
674 subnode->SetLineColor(kGreen);
675 fNodes->Add(subnode);
676 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
677 subnode->SetLineColor(kGreen);
678 fNodes->Add(subnode);
679 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
680 subnode->SetLineColor(kGreen);
681 fNodes->Add(subnode);
682 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
683 subnode->SetLineColor(kGreen);
684 fNodes->Add(subnode);
685 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
686 subnode->SetLineColor(kGreen);
687 fNodes->Add(subnode);
688 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
689 subnode->SetLineColor(kGreen);
690 fNodes->Add(subnode);
694 //Float_t pos4[3]={-171,470,0};
695 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
696 //new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
697 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
698 node->SetLineColor(kColorRICH);
700 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
701 subnode->SetLineColor(kGreen);
702 fNodes->Add(subnode);
703 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
704 subnode->SetLineColor(kGreen);
705 fNodes->Add(subnode);
706 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
707 subnode->SetLineColor(kGreen);
708 fNodes->Add(subnode);
709 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
710 subnode->SetLineColor(kGreen);
711 fNodes->Add(subnode);
712 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
713 subnode->SetLineColor(kGreen);
714 fNodes->Add(subnode);
715 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
716 subnode->SetLineColor(kGreen);
717 fNodes->Add(subnode);
722 //Float_t pos5[3]={161.3999,443.3999,-165.3};
723 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
724 //new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
725 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
726 node->SetLineColor(kColorRICH);
728 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
729 subnode->SetLineColor(kGreen);
730 fNodes->Add(subnode);
731 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
732 subnode->SetLineColor(kGreen);
733 fNodes->Add(subnode);
734 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
735 subnode->SetLineColor(kGreen);
736 fNodes->Add(subnode);
737 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
738 subnode->SetLineColor(kGreen);
739 fNodes->Add(subnode);
740 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
741 subnode->SetLineColor(kGreen);
742 fNodes->Add(subnode);
743 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
744 subnode->SetLineColor(kGreen);
745 fNodes->Add(subnode);
750 //Float_t pos6[3]={0., 471.9, -165.3,};
751 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
752 //new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
753 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
754 node->SetLineColor(kColorRICH);
755 fNodes->Add(node);node->cd();
756 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
757 subnode->SetLineColor(kGreen);
758 fNodes->Add(subnode);
759 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
760 subnode->SetLineColor(kGreen);
761 fNodes->Add(subnode);
762 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
763 subnode->SetLineColor(kGreen);
764 fNodes->Add(subnode);
765 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
766 subnode->SetLineColor(kGreen);
767 fNodes->Add(subnode);
768 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
769 subnode->SetLineColor(kGreen);
770 fNodes->Add(subnode);
771 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
772 subnode->SetLineColor(kGreen);
773 fNodes->Add(subnode);
777 //Float_t pos7[3]={-161.399,443.3999,-165.3};
778 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
779 //new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
780 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
781 node->SetLineColor(kColorRICH);
783 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
784 subnode->SetLineColor(kGreen);
785 fNodes->Add(subnode);
786 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
787 subnode->SetLineColor(kGreen);
788 fNodes->Add(subnode);
789 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
790 subnode->SetLineColor(kGreen);
791 fNodes->Add(subnode);
792 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
793 subnode->SetLineColor(kGreen);
794 fNodes->Add(subnode);
795 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
796 subnode->SetLineColor(kGreen);
797 fNodes->Add(subnode);
798 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
799 subnode->SetLineColor(kGreen);
800 fNodes->Add(subnode);
805 //___________________________________________
806 void AliRICH::CreateGeometry()
809 // Create the geometry for RICH version 1
811 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
812 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
813 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
817 <img src="picts/AliRICHv1.gif">
822 <img src="picts/AliRICHv1Tree.gif">
826 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
827 AliRICHSegmentationV0* segmentation;
828 AliRICHGeometry* geometry;
829 AliRICHChamber* iChamber;
831 iChamber = &(pRICH->Chamber(0));
832 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
833 geometry=iChamber->GetGeometryModel();
836 distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
837 geometry->SetRadiatorToPads(distance);
839 //Opaque quartz thickness
840 Float_t oqua_thickness = .5;
843 //Float_t csi_length = 160*.8 + 2.6;
844 //Float_t csi_width = 144*.84 + 2*2.6;
846 Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
847 Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
849 //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());
851 Int_t *idtmed = fIdtmed->GetArray()-999;
858 // --- Define the RICH detector
859 // External aluminium box
861 par[1] = 13; //Original Settings
866 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
870 par[1] = 13; //Original Settings
875 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
877 // Air 2 (cutting the lower part of the box)
880 par[1] = 3; //Original Settings
882 gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
884 // Air 3 (cutting the lower part of the box)
887 par[1] = 3; //Original Settings
889 gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
893 par[1] = .188; //Original Settings
898 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
902 par[1] = .025; //Original Settings
907 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
910 par[0] = geometry->GetQuartzWidth()/2;
911 par[1] = geometry->GetQuartzThickness()/2;
912 par[2] = geometry->GetQuartzLength()/2;
914 par[1] = .25; //Original Settings
916 /*par[0] = geometry->GetQuartzWidth()/2;
917 par[1] = geometry->GetQuartzThickness()/2;
918 par[2] = geometry->GetQuartzLength()/2;*/
919 //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]);
920 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
922 // Spacers (cylinders)
925 par[2] = geometry->GetFreonThickness()/2;
926 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
928 // Feet (freon slabs supports)
933 gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
936 par[0] = geometry->GetQuartzWidth()/2;
938 par[2] = geometry->GetQuartzLength()/2;
940 par[1] = .2; //Original Settings
945 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
947 // Frame of opaque quartz
948 par[0] = geometry->GetOuterFreonWidth()/2;
950 par[1] = geometry->GetFreonThickness()/2;
951 par[2] = geometry->GetOuterFreonLength()/2;
954 par[1] = .5; //Original Settings
959 gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
961 par[0] = geometry->GetInnerFreonWidth()/2;
962 par[1] = geometry->GetFreonThickness()/2;
963 par[2] = geometry->GetInnerFreonLength()/2;
964 gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
966 // Little bar of opaque quartz
968 //par[1] = geometry->GetQuartzThickness()/2;
969 //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
970 //par[2] = geometry->GetInnerFreonLength()/2;
973 par[1] = .25; //Original Settings
978 //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
981 par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
982 par[1] = geometry->GetFreonThickness()/2;
983 par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
985 par[1] = .5; //Original Settings
990 gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
992 par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
993 par[1] = geometry->GetFreonThickness()/2;
994 par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
995 gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
999 par[0] = csi_width/2;
1000 par[1] = geometry->GetGapThickness()/2;
1001 //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]);
1003 par[2] = csi_length/2;
1004 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
1008 par[0] = csi_width/2;
1009 par[1] = geometry->GetProximityGapThickness()/2;
1010 //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]);
1012 par[2] = csi_length/2;
1013 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1017 par[0] = csi_width/2;
1020 par[2] = csi_length/2;
1021 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1027 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1032 par[0] = csi_width/2;
1035 gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
1037 // Ceramic pick up (base)
1039 par[0] = csi_width/2;
1042 gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
1044 // Ceramic pick up (head)
1046 par[0] = csi_width/2;
1049 gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
1051 // Aluminium supports for methane and CsI
1054 par[0] = csi_width/2;
1055 par[1] = geometry->GetGapThickness()/2 + .25;
1056 par[2] = (68.35 - csi_length/2)/2;
1057 gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
1061 par[0] = (66.3 - csi_width/2)/2;
1062 par[1] = geometry->GetGapThickness()/2 + .25;
1063 par[2] = csi_length/2 + 68.35 - csi_length/2;
1064 gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
1066 // Aluminium supports for freon
1069 par[0] = geometry->GetQuartzWidth()/2;
1071 par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
1072 gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
1076 par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
1078 par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
1079 gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
1083 par[0] = csi_width/2;
1085 par[2] = csi_length/4 -.5025;
1086 gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
1089 // Backplane supports
1096 gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
1103 gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
1110 gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
1112 // Place holes inside backplane support
1114 gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
1115 gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
1116 gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
1117 gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
1118 gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1119 gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1120 gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1121 gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1122 gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
1123 gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
1124 gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1125 gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1126 gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1127 gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1128 gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1129 gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1133 // --- Places the detectors defined with GSVOLU
1134 // Place material inside RICH
1135 gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
1136 gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1137 gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1138 gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
1139 gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY");
1142 gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
1143 gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
1144 gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
1145 gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1146 gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
1147 gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1148 gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1149 gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1150 gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1151 gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1152 gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1153 gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
1158 gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1159 gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1160 gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
1161 gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
1165 Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
1167 gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1168 gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1169 gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1170 gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1172 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1174 //Placing of the spacers inside the freon slabs
1176 Int_t nspacers = 30;
1177 //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);
1179 //printf("Nspacers: %d", nspacers);
1181 for (i = 0; i < nspacers/3; i++) {
1182 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1183 gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1186 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1187 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1188 gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
1191 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1192 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1193 gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1196 for (i = 0; i < nspacers/3; i++) {
1197 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1198 gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1201 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1202 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1203 gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
1206 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1207 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1208 gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1212 gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
1213 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
1214 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)
1215 // printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2);
1216 gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings
1217 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)
1218 //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65)
1219 //gMC->Gspos("BARR", 2, "QUAR", geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (21.65)
1220 gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
1221 gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
1222 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1223 gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
1224 printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25);
1226 // Wire support placing
1228 gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
1229 gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
1230 gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
1232 // Backplane placing
1234 gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
1235 gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
1236 gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1237 gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1238 gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1239 gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1243 gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY");
1244 gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
1248 //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);
1250 // Place RICH inside ALICE apparatus
1254 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1255 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1256 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1257 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1258 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1259 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1260 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1262 gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
1263 gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
1264 gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
1265 gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
1266 gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
1267 gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
1268 gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");*/
1270 // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm)
1272 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
1273 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
1274 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
1275 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
1276 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
1277 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
1278 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
1280 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
1282 AliMatrix(idrotm[1000], 90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
1283 AliMatrix(idrotm[1001], 90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
1284 AliMatrix(idrotm[1002], 90., 0. , 90. , 90. , 0. , 0. );
1285 AliMatrix(idrotm[1003], 90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
1286 AliMatrix(idrotm[1004], 90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
1287 AliMatrix(idrotm[1005], 90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
1288 AliMatrix(idrotm[1006], 90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
1290 gMC->Gspos("RICH", 1, "ALIC", 0. , offset*cosphi , offset*sinphi ,idrotm[1000], "ONLY");
1291 gMC->Gspos("RICH", 2, "ALIC", (offset)*sintheta , offset*costheta , 0. ,idrotm[1001], "ONLY");
1292 gMC->Gspos("RICH", 3, "ALIC", 0. , offset , 0. ,idrotm[1002], "ONLY");
1293 gMC->Gspos("RICH", 4, "ALIC", -(offset)*sintheta, offset*costheta , 0. ,idrotm[1003], "ONLY");
1294 gMC->Gspos("RICH", 5, "ALIC", (offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1004], "ONLY");
1295 gMC->Gspos("RICH", 6, "ALIC", 0. , offset*cosphi , -offset*sinphi,idrotm[1005], "ONLY");
1296 gMC->Gspos("RICH", 7, "ALIC", -(offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1006], "ONLY");
1301 //___________________________________________
1302 void AliRICH::CreateMaterials()
1305 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1306 // ORIGIN : NICK VAN EIJNDHOVEN
1307 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1308 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
1309 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1311 Int_t isxfld = gAlice->Field()->Integ();
1312 Float_t sxmgmx = gAlice->Field()->Max();
1315 /************************************Antonnelo's Values (14-vectors)*****************************************/
1317 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,
1318 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1319 Float_t rIndexQuarz[14] = { 1.528309,1.533333,
1320 1.538243,1.544223,1.550568,1.55777,
1321 1.565463,1.574765,1.584831,1.597027,
1322 1.611858,1.6277,1.6472,1.6724 };
1323 Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1324 Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1325 Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1326 Float_t abscoFreon[14] = { 179.0987,179.0987,
1327 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
1328 //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1329 // 1e-5,1e-5,1e-5,1e-5,1e-5 };
1330 Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
1331 14.177,9.282,4.0925,1.149,.3627,.10857 };
1332 Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1333 1e-5,1e-5,1e-5,1e-5,1e-5 };
1334 Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1335 1e-4,1e-4,1e-4,1e-4 };
1336 Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1338 Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1339 1e-4,1e-4,1e-4,1e-4 };
1340 Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1341 Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
1342 .17425,.1785,.1836,.1904,.1938,.221 };
1343 Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1347 /**********************************End of Antonnelo's Values**********************************/
1349 /**********************************Values from rich_media.f (31-vectors)**********************************/
1352 //Photons energy intervals
1356 ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
1357 //printf ("Energy intervals: %e\n",ppckov[i]);
1361 //Refraction index for quarz
1362 Float_t rIndexQuarz[26];
1369 Float_t ene=ppckov[i]*1e9;
1370 Float_t a=f1/(e1*e1 - ene*ene);
1371 Float_t b=f2/(e2*e2 - ene*ene);
1372 rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
1373 //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
1376 //Refraction index for opaque quarz, methane and grid
1377 Float_t rIndexOpaqueQuarz[26];
1378 Float_t rIndexMethane[26];
1379 Float_t rIndexGrid[26];
1382 rIndexOpaqueQuarz[i]=1;
1383 rIndexMethane[i]=1.000444;
1385 //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
1388 //Absorption index for freon
1389 Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
1390 179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
1391 7.069031, 4.461292, 2.028366, 1.293013, .577267, .40746, .334964, 0., 0., 0.};
1393 //Absorption index for quarz
1394 /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1395 .906,.907,.907,.907};
1396 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,
1397 215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1398 Float_t abscoQuarz[31];
1399 for (Int_t i=0;i<31;i++)
1401 Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1402 if (Xlam <= 160) abscoQuarz[i] = 0;
1403 if (Xlam > 250) abscoQuarz[i] = 1;
1406 for (Int_t j=0;j<21;j++)
1408 //printf ("Passed\n");
1409 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1411 Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1412 Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1413 abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1417 printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1420 /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1421 39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1422 17.94531, 11.88753, 5.99128, 3.83503, 2.36661, 1.53155, 1.30582, 1.08574, .8779708,
1423 .675275, 0., 0., 0.};
1425 for (Int_t i=0;i<31;i++)
1427 abscoQuarz[i] = abscoQuarz[i]/10;
1430 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,
1431 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1432 .192, .1497, .10857};
1434 //Absorption index for methane
1435 Float_t abscoMethane[26];
1438 abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1439 //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1442 //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1443 Float_t abscoOpaqueQuarz[26];
1444 Float_t abscoCsI[26];
1445 Float_t abscoGrid[26];
1446 Float_t efficAll[26];
1447 Float_t efficGrid[26];
1450 abscoOpaqueQuarz[i]=1e-5;
1455 //printf ("All must be 1: %e, %e, %e, %e, %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1458 //Efficiency for csi
1460 Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1461 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1462 0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1463 0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1467 //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1468 //UNPOLARIZED PHOTONS
1472 efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1473 //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1476 /*******************************************End of rich_media.f***************************************/
1483 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1487 Int_t nlmatmet, nlmatqua;
1488 Float_t wmatquao[2], rIndexFreon[26];
1489 Float_t aquao[2], epsil, stmin, zquao[2];
1491 Float_t radlal, densal, tmaxfd, deemax, stemax;
1492 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1494 Int_t *idtmed = fIdtmed->GetArray()-999;
1496 // --- Photon energy (GeV)
1497 // --- Refraction indexes
1498 for (i = 0; i < 26; ++i) {
1499 rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1500 //rIndexFreon[i] = 1;
1501 //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1504 // --- Detection efficiencies (quantum efficiency for CsI)
1505 // --- Define parameters for honeycomb.
1506 // Used carbon of equivalent rad. lenght
1513 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1524 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1535 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1546 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1557 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1564 // --- Parameters to include in GSMATE related to aluminium sheet
1571 // --- Glass parameters
1573 Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
1574 Float_t zglass[5]={ 6., 14., 8., 5., 11.};
1575 Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
1576 Float_t dglass=1.74;
1579 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1580 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1581 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1582 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1583 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1584 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1585 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1586 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1587 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1588 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1589 AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1590 AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
1598 AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1599 AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1600 AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1601 AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1602 AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1603 AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1604 AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1605 AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1606 AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1607 AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1608 AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1609 AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1612 gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1613 gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1614 gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1615 gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1616 gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1617 gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1618 gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1619 gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1620 gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1621 gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1622 gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1625 //___________________________________________
1627 Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1630 //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1632 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,
1633 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,
1634 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1637 Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1638 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1639 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1642 Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1643 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1644 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1647 Int_t j=Int_t(xe*10)-49;
1648 Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1649 Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1651 //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1652 //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1654 Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1655 Float_t tanin=sinin/pdoti;
1657 Float_t c1=cn*cn-ck*ck-sinin*sinin;
1658 Float_t c2=4*cn*cn*ck*ck;
1659 Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1660 Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1662 Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1663 Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1666 //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1667 //B.J. STAGG APPLIED OPTICS, 30(1991),4113
1670 Float_t lamb=1240/ene;
1673 Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1677 Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S
1678 fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1687 //__________________________________________
1688 Float_t AliRICH::AbsoCH4(Float_t x)
1691 //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1692 Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22
1693 //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1694 Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1695 const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3
1696 const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1697 Float_t pn=kPressure/760.;
1698 Float_t tn=kTemperature/273.16;
1701 // ------- METHANE CROSS SECTION -----------------
1702 // ASTROPH. J. 214, L47 (1978)
1708 if(x>=7.75 && x<=8.1)
1710 Float_t c0=-1.655279e-1;
1711 Float_t c1=6.307392e-2;
1712 Float_t c2=-8.011441e-3;
1713 Float_t c3=3.392126e-4;
1714 sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1720 while (x<=em[j] && x>=em[j+1])
1723 Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1724 sm=(sch4[j]+a*(x-em[j]))*1e-22;
1728 Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1729 Float_t abslm=1./sm/dm;
1731 // ------- ISOBUTHANE CROSS SECTION --------------
1732 // i-C4H10 (ai) abs. length from curves in
1733 // Lu-McDonald paper for BARI RICH workshop .
1734 // -----------------------------------------------------------
1743 if(x>=7.25 && x<7.375)
1749 Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2
1750 Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1755 // ---------------------------------------------------------
1757 // transmission of O2
1759 // y= path in cm, x=energy in eV
1760 // so= cross section for UV absorption in cm2
1761 // do= O2 molecular density in cm-3
1762 // ---------------------------------------------------------
1770 so=3.392709e-13 * TMath::Exp(2.864104 *x);
1776 so=2.910039e-34 * TMath::Exp(10.3337*x);
1783 Float_t a0=-73770.76;
1784 Float_t a1=46190.69;
1785 Float_t a2=-11475.44;
1786 Float_t a3=1412.611;
1787 Float_t a4=-86.07027;
1788 Float_t a5=2.074234;
1789 so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1793 Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1798 // ---------------------------------------------------------
1800 // transmission of H2O
1802 // y= path in cm, x=energy in eV
1803 // sw= cross section for UV absorption in cm2
1804 // dw= H2O molecular density in cm-3
1805 // ---------------------------------------------------------
1809 Float_t b0=29231.65;
1810 Float_t b1=-15807.74;
1811 Float_t b2=3192.926;
1812 Float_t b3=-285.4809;
1813 Float_t b4=9.533944;
1817 Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1819 Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1825 // ---------------------------------------------------------
1827 Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1833 //___________________________________________
1834 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1842 //___________________________________________
1843 void AliRICH::MakeBranch(Option_t* option, const char *file)
1845 // Create Tree branches for the RICH.
1847 const Int_t kBufferSize = 4000;
1848 char branchname[20];
1850 AliDetector::MakeBranch(option,file);
1852 const char *cH = strstr(option,"H");
1853 const char *cD = strstr(option,"D");
1854 const char *cR = strstr(option,"R");
1855 const char *cS = strstr(option,"S");
1859 sprintf(branchname,"%sCerenkov",GetName());
1860 if (fCerenkovs && gAlice->TreeH()) {
1861 //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1862 MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1863 //branch->SetAutoDelete(kFALSE);
1865 sprintf(branchname,"%sSDigits",GetName());
1866 if (fSDigits && gAlice->TreeH()) {
1867 //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1868 MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1869 //branch->SetAutoDelete(kFALSE);
1870 //printf("Making branch %sSDigits in TreeH\n",GetName());
1875 sprintf(branchname,"%sSDigits",GetName());
1876 if (fSDigits && gAlice->TreeS()) {
1877 //TBranch* branch = MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1878 MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1879 //branch->SetAutoDelete(kFALSE);
1880 //printf("Making branch %sSDigits in TreeS\n",GetName());
1886 // one branch for digits per chamber
1890 for (i=0; i<kNCH ;i++) {
1891 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1892 if (fDchambers && gAlice->TreeD()) {
1893 //TBranch* branch = MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1894 MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1895 //branch->SetAutoDelete(kFALSE);
1896 //printf("Making Branch %sDigits%d\n",GetName(),i+1);
1903 // one branch for raw clusters per chamber
1906 //printf("Called MakeBranch for TreeR\n");
1910 for (i=0; i<kNCH ;i++) {
1911 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1912 if (fRawClusters && gAlice->TreeR()) {
1913 //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1914 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1915 //branch->SetAutoDelete(kFALSE);
1919 // one branch for rec hits per chamber
1921 for (i=0; i<kNCH ;i++) {
1922 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1923 if (fRecHits1D && gAlice->TreeR()) {
1924 //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1925 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1926 //branch->SetAutoDelete(kFALSE);
1929 for (i=0; i<kNCH ;i++) {
1930 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1931 if (fRecHits3D && gAlice->TreeR()) {
1932 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
1933 //branch->SetAutoDelete(kFALSE);
1939 //___________________________________________
1940 void AliRICH::SetTreeAddress()
1942 // Set branch address for the Hits and Digits Tree.
1943 char branchname[20];
1946 AliDetector::SetTreeAddress();
1949 TTree *treeH = gAlice->TreeH();
1950 TTree *treeD = gAlice->TreeD();
1951 TTree *treeR = gAlice->TreeR();
1952 TTree *treeS = gAlice->TreeS();
1956 branch = treeH->GetBranch("RICHCerenkov");
1957 if (branch) branch->SetAddress(&fCerenkovs);
1960 branch = treeH->GetBranch("RICHSDigits");
1963 branch->SetAddress(&fSDigits);
1964 //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits);
1971 branch = treeS->GetBranch("RICHSDigits");
1974 branch->SetAddress(&fSDigits);
1975 //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits);
1982 for (int i=0; i<kNCH; i++) {
1983 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1985 branch = treeD->GetBranch(branchname);
1986 if (branch) branch->SetAddress(&((*fDchambers)[i]));
1991 for (i=0; i<kNCH; i++) {
1992 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1994 branch = treeR->GetBranch(branchname);
1995 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1999 for (i=0; i<kNCH; i++) {
2000 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
2002 branch = treeR->GetBranch(branchname);
2003 if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
2007 for (i=0; i<kNCH; i++) {
2008 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
2010 branch = treeR->GetBranch(branchname);
2011 if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
2017 //___________________________________________
2018 void AliRICH::ResetHits()
2020 // Reset number of clusters and the cluster array for this detector
2021 AliDetector::ResetHits();
2024 if (fSDigits) fSDigits->Clear();
2025 if (fCerenkovs) fCerenkovs->Clear();
2029 //____________________________________________
2030 void AliRICH::ResetDigits()
2033 // Reset number of digits and the digits array for this detector
2035 for ( int i=0;i<kNCH;i++ ) {
2036 //PH if (fDchambers && (*fDchambers)[i]) (*fDchambers)[i]->Clear();
2037 if (fDchambers && fDchambers->At(i)) fDchambers->At(i)->Clear();
2038 if (fNdch) fNdch[i]=0;
2042 //____________________________________________
2043 void AliRICH::ResetRawClusters()
2046 // Reset number of raw clusters and the raw clust array for this detector
2048 for ( int i=0;i<kNCH;i++ ) {
2049 //PH if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
2050 if (fRawClusters->At(i)) ((TClonesArray*)fRawClusters->At(i))->Clear();
2051 if (fNrawch) fNrawch[i]=0;
2055 //____________________________________________
2056 void AliRICH::ResetRecHits1D()
2059 // Reset number of raw clusters and the raw clust array for this detector
2062 for ( int i=0;i<kNCH;i++ ) {
2063 //PH if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear();
2064 if (fRecHits1D->At(i)) ((TClonesArray*)fRecHits1D->At(i))->Clear();
2065 if (fNrechits1D) fNrechits1D[i]=0;
2069 //____________________________________________
2070 void AliRICH::ResetRecHits3D()
2073 // Reset number of raw clusters and the raw clust array for this detector
2076 for ( int i=0;i<kNCH;i++ ) {
2077 //PH if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear();
2078 if (fRecHits3D->At(i)) ((TClonesArray*)fRecHits3D->At(i))->Clear();
2079 if (fNrechits3D) fNrechits3D[i]=0;
2083 //___________________________________________
2084 void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
2088 // Setter for the RICH geometry model
2092 //PH ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
2093 ((AliRICHChamber*)fChambers->At(id))->GeometryModel(geometry);
2096 //___________________________________________
2097 void AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation)
2101 // Setter for the RICH segmentation model
2104 //PH ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation);
2105 ((AliRICHChamber*)fChambers->At(id))->SetSegmentationModel(segmentation);
2108 //___________________________________________
2109 void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
2113 // Setter for the RICH response model
2116 //PH ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
2117 ((AliRICHChamber*)fChambers->At(id))->ResponseModel(response);
2120 void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
2124 // Setter for the RICH reconstruction model (clusters)
2127 //PH ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
2128 ((AliRICHChamber*)fChambers->At(id))->SetReconstructionModel(reconst);
2131 //___________________________________________
2132 void AliRICH::StepManager()
2135 // Full Step Manager
2139 static Int_t vol[2];
2141 static Float_t hits[22];
2142 static Float_t ckovData[19];
2143 TLorentzVector position;
2144 TLorentzVector momentum;
2147 Float_t localPos[3];
2148 Float_t localMom[4];
2149 Float_t localTheta,localPhi;
2151 Float_t destep, step;
2154 Float_t coscerenkov;
2155 static Float_t eloss, xhit, yhit, tlength;
2156 const Float_t kBig=1.e10;
2158 TClonesArray &lhits = *fHits;
2159 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
2161 //if (current->Energy()>1)
2164 // Only gas gap inside chamber
2165 // Tag chambers and record hits when track enters
2168 id=gMC->CurrentVolID(copy);
2169 Float_t cherenkovLoss=0;
2170 //gAlice->KeepTrack(gAlice->CurrentTrack());
2172 gMC->TrackPosition(position);
2176 //bzero((char *)ckovData,sizeof(ckovData)*19);
2177 ckovData[1] = pos[0]; // X-position for hit
2178 ckovData[2] = pos[1]; // Y-position for hit
2179 ckovData[3] = pos[2]; // Z-position for hit
2180 ckovData[6] = 0; // dummy track length
2181 //ckovData[11] = gAlice->CurrentTrack();
2183 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
2185 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
2187 /********************Store production parameters for Cerenkov photons************************/
2188 //is it a Cerenkov photon?
2189 if (gMC->TrackPid() == 50000050) {
2191 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
2193 Float_t ckovEnergy = current->Energy();
2194 //energy interval for tracking
2195 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
2196 //if (ckovEnergy > 0)
2198 if (gMC->IsTrackEntering()){ //is track entering?
2199 //printf("Track entered (1)\n");
2200 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2202 if (gMC->IsNewTrack()){ //is it the first step?
2203 //printf("I'm in!\n");
2204 Int_t mother = current->GetFirstMother();
2206 //printf("Second Mother:%d\n",current->GetSecondMother());
2208 ckovData[10] = mother;
2209 ckovData[11] = gAlice->CurrentTrack();
2210 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
2211 //printf("Produced in FREO\n");
2214 //printf("Index: %d\n",fCkovNumber);
2215 } //first step question
2218 if (gMC->IsNewTrack()){ //is it first step?
2219 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
2222 //printf("Produced in QUAR\n");
2224 } //first step question
2226 //printf("Before %d\n",fFreonProd);
2227 } //track entering question
2229 if (ckovData[12] == 1) //was it produced in Freon?
2230 //if (fFreonProd == 1)
2232 if (gMC->IsTrackEntering()){ //is track entering?
2233 //printf("Track entered (2)\n");
2234 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
2235 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
2236 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
2238 //printf("Got in META\n");
2239 gMC->TrackMomentum(momentum);
2244 // Z-position for hit
2247 /**************** Photons lost in second grid have to be calculated by hand************/
2249 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2250 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
2252 //printf("grid calculation:%f\n",t);
2256 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2257 //printf("Added One (1)!\n");
2258 //printf("Lost one in grid\n");
2260 /**********************************************************************************/
2263 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
2264 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2265 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
2267 //printf("Got in CSI\n");
2268 gMC->TrackMomentum(momentum);
2274 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
2275 /***********************Cerenkov phtons (always polarised)*************************/
2277 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2278 Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
2283 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2284 //printf("Added One (2)!\n");
2285 //printf("Lost by Fresnel\n");
2287 /**********************************************************************************/
2292 /********************Evaluation of losses************************/
2293 /******************still in the old fashion**********************/
2296 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
2297 for (Int_t i = 0; i < i1; ++i) {
2299 if (procs[i] == kPLightReflection) { //was it reflected
2301 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2303 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2306 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2307 } //reflection question
2310 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
2311 //printf("Got in absorption\n");
2313 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2315 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2317 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
2319 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
2322 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
2326 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
2330 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2331 //printf("Added One (3)!\n");
2332 //printf("Added cerenkov %d\n",fCkovNumber);
2333 } //absorption question
2336 // Photon goes out of tracking scope
2337 else if (procs[i] == kPStop) { //is it below energy treshold?
2340 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2341 //printf("Added One (4)!\n");
2342 } // energy treshold question
2343 } //number of mechanisms cycle
2344 /**********************End of evaluation************************/
2345 } //freon production question
2346 } //energy interval question
2347 //}//inside the proximity gap question
2348 } //cerenkov photon question
2350 /**************************************End of Production Parameters Storing*********************/
2353 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
2355 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
2356 //printf("Cerenkov\n");
2358 //if (gMC->TrackPid() == 50000051)
2359 //printf("Tracking a feedback\n");
2361 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
2363 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
2364 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2365 //printf("Got in CSI\n");
2366 //printf("Tracking a %d\n",gMC->TrackPid());
2367 if (gMC->Edep() > 0.){
2368 gMC->TrackPosition(position);
2369 gMC->TrackMomentum(momentum);
2377 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2378 Double_t rt = TMath::Sqrt(tc);
2379 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2380 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2381 gMC->Gmtod(pos,localPos,1);
2382 gMC->Gmtod(mom,localMom,2);
2384 gMC->CurrentVolOffID(2,copy);
2388 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2389 //->Sector(localPos[0], localPos[2]);
2390 //printf("Sector:%d\n",sector);
2392 /*if (gMC->TrackPid() == 50000051){
2394 printf("Feedbacks:%d\n",fFeedbacks);
2397 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2398 ((AliRICHChamber*)fChambers->At(idvol))
2399 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2401 ckovData[0] = gMC->TrackPid(); // particle type
2402 ckovData[1] = pos[0]; // X-position for hit
2403 ckovData[2] = pos[1]; // Y-position for hit
2404 ckovData[3] = pos[2]; // Z-position for hit
2405 ckovData[4] = theta; // theta angle of incidence
2406 ckovData[5] = phi; // phi angle of incidence
2407 ckovData[8] = (Float_t) fNSDigits; // first sdigit
2408 ckovData[9] = -1; // last pad hit
2409 ckovData[13] = 4; // photon was detected
2410 ckovData[14] = mom[0];
2411 ckovData[15] = mom[1];
2412 ckovData[16] = mom[2];
2414 destep = gMC->Edep();
2415 gMC->SetMaxStep(kBig);
2416 cherenkovLoss += destep;
2417 ckovData[7]=cherenkovLoss;
2419 nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
2421 if (fNSDigits > (Int_t)ckovData[8]) {
2422 ckovData[8]= ckovData[8]+1;
2423 ckovData[9]= (Float_t) fNSDigits;
2426 //printf("Cerenkov loss: %f\n", cherenkovLoss);
2428 ckovData[17] = nPads;
2429 //printf("nPads:%d",nPads);
2431 //TClonesArray *Hits = RICH->Hits();
2432 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
2435 mom[0] = current->Px();
2436 mom[1] = current->Py();
2437 mom[2] = current->Pz();
2438 Float_t mipPx = mipHit->fMomX;
2439 Float_t mipPy = mipHit->fMomY;
2440 Float_t mipPz = mipHit->fMomZ;
2442 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
2443 Float_t rt = TMath::Sqrt(r);
2444 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
2445 Float_t mipRt = TMath::Sqrt(mipR);
2448 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
2454 Float_t cherenkov = TMath::ACos(coscerenkov);
2455 ckovData[18]=cherenkov;
2459 AddHit(gAlice->CurrentTrack(),vol,ckovData);
2460 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2461 //printf("Added One (5)!\n");
2468 /***********************************************End of photon hits*********************************************/
2471 /**********************************************Charged particles treatment*************************************/
2473 else if (gMC->TrackCharge())
2477 /*if (gMC->IsTrackEntering())
2479 hits[13]=20;//is track entering?
2481 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2483 gMC->TrackMomentum(momentum);
2494 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2495 // Get current particle id (ipart), track position (pos) and momentum (mom)
2497 gMC->CurrentVolOffID(3,copy);
2501 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2502 //->Sector(localPos[0], localPos[2]);
2503 //printf("Sector:%d\n",sector);
2505 gMC->TrackPosition(position);
2506 gMC->TrackMomentum(momentum);
2514 gMC->Gmtod(pos,localPos,1);
2515 gMC->Gmtod(mom,localMom,2);
2517 ipart = gMC->TrackPid();
2519 // momentum loss and steplength in last step
2520 destep = gMC->Edep();
2521 step = gMC->TrackStep();
2524 // record hits when track enters ...
2525 if( gMC->IsTrackEntering()) {
2526 // gMC->SetMaxStep(fMaxStepGas);
2527 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2528 Double_t rt = TMath::Sqrt(tc);
2529 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2530 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2533 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2534 Double_t localRt = TMath::Sqrt(localTc);
2535 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2536 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2538 hits[0] = Float_t(ipart); // particle type
2539 hits[1] = localPos[0]; // X-position for hit
2540 hits[2] = localPos[1]; // Y-position for hit
2541 hits[3] = localPos[2]; // Z-position for hit
2542 hits[4] = localTheta; // theta angle of incidence
2543 hits[5] = localPhi; // phi angle of incidence
2544 hits[8] = (Float_t) fNSDigits; // first sdigit
2545 hits[9] = -1; // last pad hit
2546 hits[13] = fFreonProd; // did id hit the freon?
2550 hits[18] = 0; // dummy cerenkov angle
2556 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2559 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2562 // Only if not trigger chamber
2565 // Initialize hit position (cursor) in the segmentation model
2566 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2567 ((AliRICHChamber*)fChambers->At(idvol))
2568 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2573 // Calculate the charge induced on a pad (disintegration) in case
2575 // Mip left chamber ...
2576 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2577 gMC->SetMaxStep(kBig);
2582 // Only if not trigger chamber
2586 if(gMC->TrackPid() == kNeutron)
2587 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2588 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2590 //printf("nPads:%d",nPads);
2596 if (fNSDigits > (Int_t)hits[8]) {
2598 hits[9]= (Float_t) fNSDigits;
2602 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2605 // Check additional signal generation conditions
2606 // defined by the segmentation
2607 // model (boundary crossing conditions)
2609 //PH (((AliRICHChamber*) (*fChambers)[idvol])
2610 (((AliRICHChamber*)fChambers->At(idvol))
2611 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2613 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2614 ((AliRICHChamber*)fChambers->At(idvol))
2615 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2618 if(gMC->TrackPid() == kNeutron)
2619 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2620 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2622 //printf("Npads:%d",NPads);
2629 // nothing special happened, add up energy loss
2636 /*************************************************End of MIP treatment**************************************/
2640 void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
2644 // Loop on chambers and on cathode planes
2646 for (Int_t icat=1;icat<2;icat++) {
2647 gAlice->ResetDigits();
2648 gAlice->TreeD()->GetEvent(0);
2649 for (Int_t ich=0;ich<kNCH;ich++) {
2650 //PH AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
2651 AliRICHChamber* iChamber=(AliRICHChamber*)fChambers->At(ich);
2652 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
2653 if (pRICHdigits == 0)
2656 // Get ready the current chamber stuff
2658 AliRICHResponse* response = iChamber->GetResponseModel();
2659 AliSegmentation* seg = iChamber->GetSegmentationModel();
2660 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2662 rec->SetSegmentation(seg);
2663 rec->SetResponse(response);
2664 rec->SetDigits(pRICHdigits);
2665 rec->SetChamber(ich);
2666 if (nev==0) rec->CalibrateCOG();
2667 rec->FindRawClusters();
2670 fRch=RawClustAddress(ich);
2674 gAlice->TreeR()->Fill();
2676 for (int i=0;i<kNCH;i++) {
2677 fRch=RawClustAddress(i);
2678 int nraw=fRch->GetEntriesFast();
2679 printf ("Chamber %d, raw clusters %d\n",i,nraw);
2687 sprintf(hname,"TreeR%d",nev);
2688 gAlice->TreeR()->Write(hname,kOverwrite,0);
2689 gAlice->TreeR()->Reset();
2691 //gObjectTable->Print();
2694 AliRICHSDigit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
2697 // Initialise the pad iterator
2698 // Return the address of the first sdigit for hit
2699 TClonesArray *theClusters = clusters;
2700 Int_t nclust = theClusters->GetEntriesFast();
2701 if (nclust && hit->fPHlast > 0) {
2702 sMaxIterPad=Int_t(hit->fPHlast);
2703 sCurIterPad=Int_t(hit->fPHfirst);
2704 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2711 AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters)
2714 // Iterates over pads
2717 if (sCurIterPad <= sMaxIterPad) {
2718 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2724 AliRICH& AliRICH::operator=(const AliRICH& rhs)
2726 // Assignment operator
2731 void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
2734 Int_t NpadX = 162; // number of pads on X
2735 Int_t NpadY = 162; // number of pads on Y
2737 Int_t Pad[162][162];
2738 for (Int_t i=0;i<NpadX;i++) {
2739 for (Int_t j=0;j<NpadY;j++) {
2744 // Create some histograms
2746 TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
2747 TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
2748 TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
2749 TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
2750 TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
2751 TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
2752 TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
2753 TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
2754 TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
2755 TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
2756 TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
2757 TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
2758 TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
2759 TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
2760 TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
2761 TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
2762 TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2763 TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2764 TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
2765 TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2766 TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2767 TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
2768 TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
2769 TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
2770 TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
2771 //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
2772 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
2773 TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
2774 TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
2775 TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
2776 TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
2777 TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
2782 // Start loop over events
2784 Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
2785 Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
2786 Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
2789 for (int nev=0; nev<= evNumber2; nev++) {
2790 Int_t nparticles = gAlice->GetEvent(nev);
2793 printf ("Event number : %d\n",nev);
2794 printf ("Number of particles: %d\n",nparticles);
2795 if (nev < evNumber1) continue;
2796 if (nparticles <= 0) return;
2798 // Get pointers to RICH detector and Hits containers
2800 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
2802 TTree *treeH = gAlice->TreeH();
2803 Int_t ntracks =(Int_t) treeH->GetEntries();
2805 // Start loop on tracks in the hits containers
2807 for (Int_t track=0; track<ntracks;track++) {
2808 printf ("Processing Track: %d\n",track);
2809 gAlice->ResetHits();
2810 treeH->GetEvent(track);
2812 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2814 mHit=(AliRICHHit*)pRICH->NextHit())
2816 //Int_t nch = mHit->fChamber; // chamber number
2817 //Float_t x = mHit->X(); // x-pos of hit
2818 //Float_t y = mHit->Z(); // y-pos
2819 //Float_t z = mHit->Y();
2820 //Float_t phi = mHit->fPhi; //Phi angle of incidence
2821 Float_t theta = mHit->fTheta; //Theta angle of incidence
2822 Float_t px = mHit->MomX();
2823 Float_t py = mHit->MomY();
2824 Int_t index = mHit->Track();
2825 Int_t particle = (Int_t)(mHit->fParticle);
2830 TParticle *current = gAlice->Particle(index);
2832 //Float_t energy=current->Energy();
2834 R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
2835 PTfinal=TMath::Sqrt(px*px + py*py);
2836 PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
2840 if (TMath::Abs(particle) < 10000000)
2842 hitsTheta->Fill(theta,(float) 1);
2845 if (PTvertex>.5 && PTvertex<=1)
2847 hitsTheta500MeV->Fill(theta,(float) 1);
2849 if (PTvertex>1 && PTvertex<=2)
2851 hitsTheta1GeV->Fill(theta,(float) 1);
2853 if (PTvertex>2 && PTvertex<=3)
2855 hitsTheta2GeV->Fill(theta,(float) 1);
2859 hitsTheta3GeV->Fill(theta,(float) 1);
2868 //printf("Particle type: %d\n",current->GetPdgCode());
2869 if (TMath::Abs(particle) < 50000051)
2871 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
2872 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
2874 //gMC->Rndm(&random, 1);
2875 if (random->Rndm() < .1)
2876 production->Fill(current->Vz(),R,(float) 1);
2877 if (TMath::Abs(particle) == 50000050)
2878 //if (TMath::Abs(particle) > 50000000)
2884 if (current->Energy()>0.001)
2885 highprimaryphotons +=1;
2888 if (TMath::Abs(particle) == 2112)
2891 if (current->Energy()>0.0001)
2895 if (TMath::Abs(particle) < 50000000)
2897 production->Fill(current->Vz(),R,(float) 1);
2898 //printf("Adding %d at %f\n",particle,R);
2900 //mip->Fill(x,y,(float) 1);
2903 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2907 pionptspectravertex->Fill(PTvertex,(float) 1);
2908 pionptspectrafinal->Fill(PTfinal,(float) 1);
2912 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2913 || TMath::Abs(particle)==311)
2917 kaonptspectravertex->Fill(PTvertex,(float) 1);
2918 kaonptspectrafinal->Fill(PTfinal,(float) 1);
2923 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2925 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2926 //printf ("fParticle: %d, PDG code:%d\n",particle,current->GetPdgCode());
2927 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2928 pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2931 pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2932 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2934 //printf("Pion mass: %e\n",current->GetCalcMass());
2936 if (TMath::Abs(particle)==211)
2942 if (current->Energy()>1)
2943 highprimarypions +=1;
2947 if (TMath::Abs(particle)==2212)
2949 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2950 //ptspectra->Fill(Pt,(float) 1);
2951 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2952 protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2954 protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2955 //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
2958 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2959 || TMath::Abs(particle)==311)
2961 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2962 //ptspectra->Fill(Pt,(float) 1);
2963 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2964 kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2966 kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2967 //printf("Kaon mass: %e\n",current->GetCalcMass());
2969 if (TMath::Abs(particle)==321)
2975 if (current->Energy()>1)
2976 highprimarykaons +=1;
2980 if (TMath::Abs(particle)==11)
2982 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2983 //ptspectra->Fill(Pt,(float) 1);
2984 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2985 electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2987 electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2988 //printf("Electron mass: %e\n",current->GetCalcMass());
2991 if (particle == -11)
2994 if (TMath::Abs(particle)==13)
2996 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2997 //ptspectra->Fill(Pt,(float) 1);
2998 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2999 muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3001 muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3002 //printf("Muon mass: %e\n",current->GetCalcMass());
3005 if (TMath::Abs(particle)==2112)
3007 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3008 //ptspectra->Fill(Pt,(float) 1);
3009 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
3010 neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3013 neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3014 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
3016 //printf("Neutron mass: %e\n",current->GetCalcMass());
3019 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3021 if (current->Energy()-current->GetCalcMass()>1)
3023 chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3024 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
3025 chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3027 chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3030 //printf("Hits:%d\n",hit);
3031 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3032 // Fill the histograms
3034 //h->Fill(x,y,(float) 1);
3044 //Create canvases, set the view range, show histograms
3046 TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
3048 //c2->SetFillColor(42);
3051 hitsTheta500MeV->SetFillColor(5);
3052 hitsTheta500MeV->Draw();
3054 hitsTheta1GeV->SetFillColor(5);
3055 hitsTheta1GeV->Draw();
3057 hitsTheta2GeV->SetFillColor(5);
3058 hitsTheta2GeV->Draw();
3060 hitsTheta3GeV->SetFillColor(5);
3061 hitsTheta3GeV->Draw();
3065 TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
3067 production->SetFillColor(42);
3068 production->SetXTitle("z (m)");
3069 production->SetYTitle("R (m)");
3072 TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
3075 pionptspectravertex->SetFillColor(5);
3076 pionptspectravertex->SetXTitle("Pt (GeV)");
3077 pionptspectravertex->Draw();
3079 pionptspectrafinal->SetFillColor(5);
3080 pionptspectrafinal->SetXTitle("Pt (GeV)");
3081 pionptspectrafinal->Draw();
3083 kaonptspectravertex->SetFillColor(5);
3084 kaonptspectravertex->SetXTitle("Pt (GeV)");
3085 kaonptspectravertex->Draw();
3087 kaonptspectrafinal->SetFillColor(5);
3088 kaonptspectrafinal->SetXTitle("Pt (GeV)");
3089 kaonptspectrafinal->Draw();
3092 TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
3096 //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
3097 electronspectra1->SetFillColor(5);
3098 electronspectra1->SetXTitle("log(GeV)");
3099 electronspectra2->SetFillColor(46);
3100 electronspectra2->SetXTitle("log(GeV)");
3101 electronspectra3->SetFillColor(10);
3102 electronspectra3->SetXTitle("log(GeV)");
3104 electronspectra1->Draw();
3105 electronspectra2->Draw("same");
3106 electronspectra3->Draw("same");
3109 //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
3110 muonspectra1->SetFillColor(5);
3111 muonspectra1->SetXTitle("log(GeV)");
3112 muonspectra2->SetFillColor(46);
3113 muonspectra2->SetXTitle("log(GeV)");
3114 muonspectra3->SetFillColor(10);
3115 muonspectra3->SetXTitle("log(GeV)");
3117 muonspectra1->Draw();
3118 muonspectra2->Draw("same");
3119 muonspectra3->Draw("same");
3122 //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
3123 //neutronspectra1->SetFillColor(42);
3124 //neutronspectra1->SetXTitle("log(GeV)");
3125 //neutronspectra2->SetFillColor(46);
3126 //neutronspectra2->SetXTitle("log(GeV)");
3127 //neutronspectra3->SetFillColor(10);
3128 //neutronspectra3->SetXTitle("log(GeV)");
3130 //neutronspectra1->Draw();
3131 //neutronspectra2->Draw("same");
3132 //neutronspectra3->Draw("same");
3134 TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
3135 //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
3139 pionspectra1->SetFillColor(5);
3140 pionspectra1->SetXTitle("log(GeV)");
3141 pionspectra2->SetFillColor(46);
3142 pionspectra2->SetXTitle("log(GeV)");
3143 pionspectra3->SetFillColor(10);
3144 pionspectra3->SetXTitle("log(GeV)");
3146 pionspectra1->Draw();
3147 pionspectra2->Draw("same");
3148 pionspectra3->Draw("same");
3151 //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
3152 protonspectra1->SetFillColor(5);
3153 protonspectra1->SetXTitle("log(GeV)");
3154 protonspectra2->SetFillColor(46);
3155 protonspectra2->SetXTitle("log(GeV)");
3156 protonspectra3->SetFillColor(10);
3157 protonspectra3->SetXTitle("log(GeV)");
3159 protonspectra1->Draw();
3160 protonspectra2->Draw("same");
3161 protonspectra3->Draw("same");
3164 //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
3165 kaonspectra1->SetFillColor(5);
3166 kaonspectra1->SetXTitle("log(GeV)");
3167 kaonspectra2->SetFillColor(46);
3168 kaonspectra2->SetXTitle("log(GeV)");
3169 kaonspectra3->SetFillColor(10);
3170 kaonspectra3->SetXTitle("log(GeV)");
3172 kaonspectra1->Draw();
3173 kaonspectra2->Draw("same");
3174 kaonspectra3->Draw("same");
3177 //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
3178 chargedspectra1->SetFillColor(5);
3179 chargedspectra1->SetXTitle("log(GeV)");
3180 chargedspectra2->SetFillColor(46);
3181 chargedspectra2->SetXTitle("log(GeV)");
3182 chargedspectra3->SetFillColor(10);
3183 chargedspectra3->SetXTitle("log(GeV)");
3185 chargedspectra1->Draw();
3186 chargedspectra2->Draw("same");
3187 chargedspectra3->Draw("same");
3191 printf("*****************************************\n");
3192 printf("* Particle * Counts *\n");
3193 printf("*****************************************\n");
3195 printf("* Pions: * %4d *\n",pion);
3196 printf("* Charged Pions: * %4d *\n",chargedpions);
3197 printf("* Primary Pions: * %4d *\n",primarypions);
3198 printf("* Primary Pions (p>1GeV/c): * %4d *\n",highprimarypions);
3199 printf("* Kaons: * %4d *\n",kaon);
3200 printf("* Charged Kaons: * %4d *\n",chargedkaons);
3201 printf("* Primary Kaons: * %4d *\n",primarykaons);
3202 printf("* Primary Kaons (p>1GeV/c): * %4d *\n",highprimarykaons);
3203 printf("* Muons: * %4d *\n",muon);
3204 printf("* Electrons: * %4d *\n",electron);
3205 printf("* Positrons: * %4d *\n",positron);
3206 printf("* Protons: * %4d *\n",proton);
3207 printf("* All Charged: * %4d *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
3208 printf("*****************************************\n");
3209 //printf("* Photons: * %3.1f *\n",photons);
3210 //printf("* Primary Photons: * %3.1f *\n",primaryphotons);
3211 //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons);
3212 //printf("*****************************************\n");
3213 //printf("* Neutrons: * %3.1f *\n",neutron);
3214 //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons);
3215 //printf("*****************************************\n");
3217 if (gAlice->TreeD())
3219 gAlice->TreeD()->GetEvent(0);
3224 printf("\n*****************************************\n");
3225 printf("* Chamber * Digits * Occupancy *\n");
3226 printf("*****************************************\n");
3228 for (Int_t ich=0;ich<7;ich++)
3230 TClonesArray *Digits = DigitsAddress(ich); // Raw clusters branch
3231 Int_t ndigits = Digits->GetEntriesFast();
3232 occ[ich] = Float_t(ndigits)/(160*144);
3233 sum += Float_t(ndigits)/(160*144);
3234 printf("* %d * %d * %3.1f%% *\n",ich,ndigits,occ[ich]*100);
3237 printf("*****************************************\n");
3238 printf("* Mean occupancy * %3.1f%% *\n",mean*100);
3239 printf("*****************************************\n");
3242 printf("\nEnd of analysis\n");
3246 //_________________________________________________________________________________________________
3249 void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
3252 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
3253 AliRICHSegmentationV0* segmentation;
3254 AliRICHChamber* chamber;
3256 chamber = &(pRICH->Chamber(0));
3257 segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel(0);
3259 Int_t NpadX = segmentation->Npx(); // number of pads on X
3260 Int_t NpadY = segmentation->Npy(); // number of pads on Y
3262 //Int_t Pad[144][160];
3263 /*for (Int_t i=0;i<NpadX;i++) {
3264 for (Int_t j=0;j<NpadY;j++) {
3270 Int_t xmin= -NpadX/2;
3271 Int_t xmax= NpadX/2;
3272 Int_t ymin= -NpadY/2;
3273 Int_t ymax= NpadY/2;
3282 TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-30,30,150,-50,10);
3286 printf("Single Ring Hits\n");
3287 feedback = new TH2F("feedback","Feedback hit distribution",150,-30,30,150,-50,10);
3288 mip = new TH2F("mip","Mip hit distribution",150,-30,30,150,-50,10);
3289 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-30,30,150,-50,10);
3290 h = new TH2F("h","Detector hit distribution",150,-30,30,150,-50,10);
3291 hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-30,30);
3292 hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,10);
3296 printf("Full Event Hits\n");
3298 feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
3299 mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
3300 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
3301 h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
3302 hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
3303 hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
3308 TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3309 TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3310 TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3311 TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3312 TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3313 TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3314 TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3316 TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
3317 TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.3,1);
3318 TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
3319 TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
3320 TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
3321 TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
3322 TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
3323 TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
3324 TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
3325 //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
3326 TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
3327 TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
3328 TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
3329 TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
3330 TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
3331 TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
3332 TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
3333 TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
3334 TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
3335 TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
3336 TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
3337 TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,0,360);
3338 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15);
3339 TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1);
3340 TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15);
3341 TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,0,360);
3342 TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.3,1);
3343 TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",200,.3,1);
3344 TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
3345 TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
3347 // Start loop over events
3352 Int_t mothers[80000];
3353 Int_t mothers2[80000];
3361 for (Int_t i=0;i<100;i++) mothers[i]=0;
3363 for (int nev=0; nev<= evNumber2; nev++) {
3364 Int_t nparticles = gAlice->GetEvent(nev);
3367 //cout<<"nev "<<nev<<endl;
3368 printf ("\n**********************************\nProcessing Event: %d\n",nev);
3369 //cout<<"nparticles "<<nparticles<<endl;
3370 printf ("Particles : %d\n\n",nparticles);
3371 if (nev < evNumber1) continue;
3372 if (nparticles <= 0) return;
3374 // Get pointers to RICH detector and Hits containers
3377 TTree *TH = gAlice->TreeH();
3378 Stat_t ntracks = TH->GetEntries();
3380 // Start loop on tracks in the hits containers
3382 for (Int_t track=0; track<ntracks;track++) {
3384 printf ("\nProcessing Track: %d\n",track);
3385 gAlice->ResetHits();
3386 TH->GetEvent(track);
3387 Int_t nhits = pRICH->Hits()->GetEntriesFast();
3388 if (nhits) Nh+=nhits;
3389 printf("Hits : %d\n",nhits);
3390 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
3392 mHit=(AliRICHHit*)pRICH->NextHit())
3394 //Int_t nch = mHit->fChamber; // chamber number
3395 x = mHit->X(); // x-pos of hit
3396 y = mHit->Z(); // y-pos
3397 Float_t phi = mHit->fPhi; //Phi angle of incidence
3398 Float_t theta = mHit->fTheta; //Theta angle of incidence
3399 Int_t index = mHit->Track();
3400 Int_t particle = (Int_t)(mHit->fParticle);
3401 //Int_t freon = (Int_t)(mHit->fLoss);
3403 hitsX->Fill(x,(float) 1);
3404 hitsY->Fill(y,(float) 1);
3406 //printf("Particle:%9d\n",particle);
3408 TParticle *current = (TParticle*)gAlice->Particle(index);
3409 //printf("Particle type: %d\n",sizeoff(Particles));
3411 hitsTheta->Fill(theta,(float) 1);
3412 //hitsPhi->Fill(phi,(float) 1);
3413 //if (pRICH->GetDebugLevel() == -1)
3414 //printf("Theta:%f, Phi:%f\n",theta,phi);
3416 //printf("Debug Level:%d\n",pRICH->GetDebugLevel());
3418 if (current->GetPdgCode() < 10000000)
3420 mip->Fill(x,y,(float) 1);
3421 //printf("adding mip\n");
3422 //if (current->Energy() - current->GetCalcMass()>1 && freon==1)
3424 hitsPhi->Fill(TMath::Abs(phi),(float) 1);
3425 //hitsTheta->Fill(theta,(float) 1);
3426 //printf("Theta:%f, Phi:%f\n",theta,phi);
3430 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
3432 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3434 if (TMath::Abs(particle)==2212)
3436 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3438 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
3439 || TMath::Abs(particle)==311)
3441 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3443 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3445 if (current->Energy() - current->GetCalcMass()>1)
3446 chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3448 //printf("Hits:%d\n",hit);
3449 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3450 // Fill the histograms
3452 h->Fill(x,y,(float) 1);
3457 Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
3458 //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
3459 //totalphotonsevent->Fill(ncerenkovs,(float) 1);
3462 printf("Cerenkovs : %d\n",ncerenkovs);
3463 totalphotonsevent->Fill(ncerenkovs,(float) 1);
3464 for (Int_t hit=0;hit<ncerenkovs;hit++) {
3465 AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
3466 //Int_t nchamber = cHit->fChamber; // chamber number
3467 Int_t index = cHit->Track();
3468 //Int_t pindex = (Int_t)(cHit->fIndex);
3469 Float_t cx = cHit->X(); // x-position
3470 Float_t cy = cHit->Z(); // y-position
3471 Int_t cmother = cHit->fCMother; // Index of mother particle
3472 Int_t closs = (Int_t)(cHit->fLoss); // How did the particle get lost?
3473 Float_t cherenkov = cHit->fCerenkovAngle; //production cerenkov angle
3474 //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
3476 //printf("Particle:%9d\n",index);
3478 TParticle *current = (TParticle*)gAlice->Particle(index);
3479 Float_t energyckov = current->Energy();
3481 if (current->GetPdgCode() == 50000051)
3485 feedback->Fill(cx,cy,(float) 1);
3489 if (current->GetPdgCode() == 50000050)
3494 phspectra2->Fill(energyckov*1e9,(float) 1);
3499 cerenkov->Fill(cx,cy,(float) 1);
3501 //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
3503 //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
3504 AliRICHHit* mipHit = (AliRICHHit*) pRICH->Hits()->UncheckedAt(0);
3505 mom[0] = current->Px();
3506 mom[1] = current->Py();
3507 mom[2] = current->Pz();
3508 //mom[0] = cHit->fMomX;
3509 // mom[1] = cHit->fMomZ;
3510 //mom[2] = cHit->fMomY;
3511 //Float_t energymip = MIP->Energy();
3512 //Float_t Mip_px = mipHit->fMomFreoX;
3513 //Float_t Mip_py = mipHit->fMomFreoY;
3514 //Float_t Mip_pz = mipHit->fMomFreoZ;
3515 //Float_t Mip_px = MIP->Px();
3516 //Float_t Mip_py = MIP->Py();
3517 //Float_t Mip_pz = MIP->Pz();
3521 //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
3522 //Float_t rt = TMath::Sqrt(r);
3523 //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
3524 //Float_t Mip_rt = TMath::Sqrt(Mip_r);
3525 //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
3526 //Float_t cherenkov = TMath::ACos(coscerenkov);
3527 ckovangle->Fill(cherenkov,(float) 1); //Cerenkov angle calculus
3528 //printf("Cherenkov: %f\n",cherenkov);
3529 Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
3530 hckphi->Fill(ckphi,(float) 1);
3533 //Float_t mix = MIP->Vx();
3534 //Float_t miy = MIP->Vy();
3535 Float_t mx = mipHit->X();
3536 Float_t my = mipHit->Z();
3537 //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
3538 Float_t dx = cx - mx;
3539 Float_t dy = cy - my;
3540 //printf("Dx:%f, Dy:%f\n",dx,dy);
3541 Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
3542 //printf("Final radius:%f\n",final_radius);
3543 radius->Fill(final_radius,(float) 1);
3545 phspectra1->Fill(energyckov*1e9,(float) 1);
3548 for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
3549 if (cmother == nmothers){
3551 mothers2[cmother]++;
3562 Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
3563 gAlice->TreeR()->GetEvent(nent-1);
3564 TClonesArray *Rawclusters = pRICH->RawClustAddress(2); // Raw clusters branch
3565 //printf ("Rawclusters:%p",Rawclusters);
3566 Int_t nrawclusters = Rawclusters->GetEntriesFast();
3569 printf("Raw Clusters : %d\n",nrawclusters);
3570 for (Int_t hit=0;hit<nrawclusters;hit++) {
3571 AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
3572 //Int_t nchamber = rcHit->fChamber; // chamber number
3573 //Int_t nhit = cHit->fHitNumber; // hit number
3574 Int_t qtot = rcHit->fQ; // charge
3575 Float_t fx = rcHit->fX; // x-position
3576 Float_t fy = rcHit->fY; // y-position
3577 //Int_t type = rcHit->fCtype; // cluster type ?
3578 Int_t mult = rcHit->fMultiplicity; // How many pads form the cluster
3581 //printf ("fx: %d, fy: %d\n",fx,fy);
3582 if (fx>(x-4) && fx<(x+4) && fy>(y-4) && fy<(y+4)) {
3583 //printf("There %d \n",mult);
3586 padnumber->Fill(mult,(float) 1);
3588 if (mult<4) Clcharge->Fill(qtot,(float) 1);
3596 TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3597 Int_t nrechits1D = RecHits1D->GetEntriesFast();
3598 //printf (" nrechits:%d\n",nrechits);
3602 for (Int_t hit=0;hit<nrechits1D;hit++) {
3603 AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
3604 Float_t r_omega = recHit1D->fOmega; // Cerenkov angle
3605 Float_t *cer_pho = recHit1D->fCerPerPhoton; // Cerenkov angle per photon
3606 Int_t *padsx = recHit1D->fPadsUsedX; // Pads Used fo reconstruction (x)
3607 Int_t *padsy = recHit1D->fPadsUsedY; // Pads Used fo reconstruction (y)
3608 Int_t goodPhotons = recHit1D->fGoodPhotons; // Number of pads used for reconstruction
3610 Omega1D->Fill(r_omega,(float) 1);
3612 for (Int_t i=0; i<goodPhotons; i++)
3614 PhotonCer->Fill(cer_pho[i],(float) 1);
3615 PadsUsed->Fill(padsx[i],padsy[i],1);
3616 //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
3619 //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
3624 TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3625 Int_t nrechits3D = RecHits3D->GetEntriesFast();
3626 //printf (" nrechits:%d\n",nrechits);
3630 for (Int_t hit=0;hit<nrechits3D;hit++) {
3631 AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(hit);
3632 Float_t r_omega = recHit3D->fOmega; // Cerenkov angle
3633 Float_t r_theta = recHit3D->fTheta; // Theta angle of incidence
3634 Float_t r_phi = recHit3D->fPhi; // Phi angle if incidence
3635 Float_t meanradius = recHit3D->fMeanRadius; // Mean radius for reconstructed point
3637 //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY);
3639 Omega3D->Fill(r_omega,(float) 1);
3640 Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
3641 Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
3642 MeanRadius->Fill(meanradius,(float) 1);
3648 for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
3649 totalphotonstrack->Fill(mothers[nmothers],(float) 1);
3650 mother->Fill(mothers2[nmothers],(float) 1);
3651 //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
3654 clusev->Fill(nraw,(float) 1);
3655 photev->Fill(phot,(float) 1);
3656 feedev->Fill(feed,(float) 1);
3657 padsmip->Fill(padmip,(float) 1);
3658 padscl->Fill(pads,(float) 1);
3659 //printf("Photons:%d\n",phot);
3668 gAlice->ResetDigits();
3669 //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
3670 gAlice->TreeD()->GetEvent(0);
3676 TClonesArray *Digits = pRICH->DigitsAddress(2);
3677 Int_t ndigits = Digits->GetEntriesFast();
3678 printf("Digits : %d\n",ndigits);
3679 padsev->Fill(ndigits,(float) 1);
3680 for (Int_t hit=0;hit<ndigits;hit++) {
3681 AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3682 Int_t qtot = dHit->fSignal; // charge
3683 Int_t ipx = dHit->fPadX; // pad number on X
3684 Int_t ipy = dHit->fPadY; // pad number on Y
3685 //printf("%d, %d\n",ipx,ipy);
3686 if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
3692 for (Int_t ich=0;ich<7;ich++)
3694 TClonesArray *Digits = pRICH->DigitsAddress(ich); // Raw clusters branch
3695 Int_t ndigits = Digits->GetEntriesFast();
3696 //printf("Digits:%d\n",ndigits);
3697 padsev->Fill(ndigits,(float) 1);
3699 for (Int_t hit=0;hit<ndigits;hit++) {
3700 AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3701 //Int_t nchamber = dHit->fChamber; // chamber number
3702 //Int_t nhit = dHit->fHitNumber; // hit number
3703 Int_t qtot = dHit->fSignal; // charge
3704 Int_t ipx = dHit->fPadX; // pad number on X
3705 Int_t ipy = dHit->fPadY; // pad number on Y
3706 //Int_t iqpad = dHit->fQpad; // charge per pad
3707 //Int_t rpad = dHit->fRSec; // R-position of pad
3708 //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
3709 if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
3710 if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
3711 if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
3712 if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
3713 if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
3714 if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
3715 if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
3716 if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
3724 //Create canvases, set the view range, show histograms
3742 TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3743 Int_t nrechits3D = RecHits3D->GetEntriesFast();
3744 TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3745 Int_t nrechits1D = RecHits1D->GetEntriesFast();
3751 c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
3752 hc0->SetXTitle("ix (npads)");
3756 c2 = new TCanvas("c2","Hits per type",100,100,600,700);
3758 //c4->SetFillColor(42);
3761 feedback->SetXTitle("x (cm)");
3762 feedback->SetYTitle("y (cm)");
3766 //mip->SetFillColor(5);
3767 mip->SetXTitle("x (cm)");
3768 mip->SetYTitle("y (cm)");
3772 //cerenkov->SetFillColor(5);
3773 cerenkov->SetXTitle("x (cm)");
3774 cerenkov->SetYTitle("y (cm)");
3778 //h->SetFillColor(5);
3779 h->SetXTitle("x (cm)");
3780 h->SetYTitle("y (cm)");
3783 c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
3785 //c10->SetFillColor(42);
3788 hitsX->SetFillColor(5);
3789 hitsX->SetXTitle("(cm)");
3793 hitsY->SetFillColor(5);
3794 hitsY->SetXTitle("(cm)");
3802 c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
3804 //c6->SetFillColor(42);
3807 phspectra2->SetFillColor(5);
3808 phspectra2->SetXTitle("energy (eV)");
3811 phspectra1->SetFillColor(5);
3812 phspectra1->SetXTitle("energy (eV)");
3815 c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
3817 //c9->SetFillColor(42);
3820 pionspectra->SetFillColor(5);
3821 pionspectra->SetXTitle("(GeV)");
3822 pionspectra->Draw();
3825 protonspectra->SetFillColor(5);
3826 protonspectra->SetXTitle("(GeV)");
3827 protonspectra->Draw();
3830 kaonspectra->SetFillColor(5);
3831 kaonspectra->SetXTitle("(GeV)");
3832 kaonspectra->Draw();
3835 chargedspectra->SetFillColor(5);
3836 chargedspectra->SetXTitle("(GeV)");
3837 chargedspectra->Draw();
3846 c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
3848 //c3->SetFillColor(42);
3853 Clcharge->SetFillColor(5);
3854 Clcharge->SetXTitle("ADC counts");
3857 Clcharge->Fit("expo");
3858 //expo->SetLineColor(2);
3859 //expo->SetLineWidth(3);
3864 padnumber->SetFillColor(5);
3865 padnumber->SetXTitle("(counts)");
3869 clusev->SetFillColor(5);
3870 clusev->SetXTitle("(counts)");
3873 clusev->Fit("gaus");
3874 //gaus->SetLineColor(2);
3875 //gaus->SetLineWidth(3);
3880 padsmip->SetFillColor(5);
3881 padsmip->SetXTitle("(counts)");
3887 c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
3888 mother->SetFillColor(5);
3889 mother->SetXTitle("counts");
3893 c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
3895 //c7->SetFillColor(42);
3898 totalphotonsevent->SetFillColor(5);
3899 totalphotonsevent->SetXTitle("Photons (counts)");
3902 totalphotonsevent->Fit("gaus");
3903 //gaus->SetLineColor(2);
3904 //gaus->SetLineWidth(3);
3906 totalphotonsevent->Draw();
3909 photev->SetFillColor(5);
3910 photev->SetXTitle("(counts)");
3913 photev->Fit("gaus");
3914 //gaus->SetLineColor(2);
3915 //gaus->SetLineWidth(3);
3920 feedev->SetFillColor(5);
3921 feedev->SetXTitle("(counts)");
3924 feedev->Fit("gaus");
3925 //gaus->SetLineColor(2);
3926 //gaus->SetLineWidth(3);
3931 padsev->SetFillColor(5);
3932 padsev->SetXTitle("(counts)");
3935 padsev->Fit("gaus");
3936 //gaus->SetLineColor(2);
3937 //gaus->SetLineWidth(3);
3948 c8 = new TCanvas("c8","3D reconstruction",50,50,1100,700);
3950 //c2->SetFillColor(42);
3953 hitsPhi->SetFillColor(5);
3956 hitsTheta->SetFillColor(5);
3959 ckovangle->SetFillColor(5);
3960 ckovangle->SetXTitle("angle (radians)");
3963 radius->SetFillColor(5);
3964 radius->SetXTitle("radius (cm)");
3967 Phi->SetFillColor(5);
3970 Theta->SetFillColor(5);
3973 Omega3D->SetFillColor(5);
3974 Omega3D->SetXTitle("angle (radians)");
3977 MeanRadius->SetFillColor(5);
3978 MeanRadius->SetXTitle("radius (cm)");
3985 c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
3987 //c5->SetFillColor(42);
3990 ckovangle->SetFillColor(5);
3991 ckovangle->SetXTitle("angle (radians)");
3995 radius->SetFillColor(5);
3996 radius->SetXTitle("radius (cm)");
4000 hc0->SetXTitle("pads");
4004 Omega1D->SetFillColor(5);
4005 Omega1D->SetXTitle("angle (radians)");
4009 PhotonCer->SetFillColor(5);
4010 PhotonCer->SetXTitle("angle (radians)");
4014 PadsUsed->SetXTitle("pads");
4015 PadsUsed->Draw("box");
4022 printf("Drawing histograms.../n");
4026 c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
4028 //c1->SetFillColor(42);
4031 hc1->SetXTitle("ix (npads)");
4034 hc2->SetXTitle("ix (npads)");
4037 hc3->SetXTitle("ix (npads)");
4040 hc4->SetXTitle("ix (npads)");
4043 hc5->SetXTitle("ix (npads)");
4046 hc6->SetXTitle("ix (npads)");
4049 hc7->SetXTitle("ix (npads)");
4052 hc0->SetXTitle("ix (npads)");
4056 c11 = new TCanvas("c11","Hits per type",100,100,600,700);
4058 //c4->SetFillColor(42);
4061 feedback->SetXTitle("x (cm)");
4062 feedback->SetYTitle("y (cm)");
4066 //mip->SetFillColor(5);
4067 mip->SetXTitle("x (cm)");
4068 mip->SetYTitle("y (cm)");
4072 //cerenkov->SetFillColor(5);
4073 cerenkov->SetXTitle("x (cm)");
4074 cerenkov->SetYTitle("y (cm)");
4078 //h->SetFillColor(5);
4079 h->SetXTitle("x (cm)");
4080 h->SetYTitle("y (cm)");
4083 c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
4085 //c10->SetFillColor(42);
4088 hitsX->SetFillColor(5);
4089 hitsX->SetXTitle("(cm)");
4093 hitsY->SetFillColor(5);
4094 hitsY->SetXTitle("(cm)");
4102 // calculate the number of pads which give a signal
4106 /*for (Int_t i=0;i< NpadX;i++) {
4107 for (Int_t j=0;j< NpadY;j++) {
4113 //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
4114 printf("\nEnd of analysis\n");
4115 printf("**********************************\n");