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