New full version. All parameters configurable.
[u/mrichter/AliRoot.git] / RICH / AliRICHv1.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 /*
17   $Log$
18   Revision 1.9  2000/05/31 08:19:38  jbarbosa
19   Fixed bug in StepManager
20
21   Revision 1.8  2000/05/26 17:30:08  jbarbosa
22   Cerenkov angle now stored within cerenkov data structure.
23
24   Revision 1.7  2000/05/18 10:31:36  jbarbosa
25   Fixed positioning of spacers inside freon.
26   Fixed positioning of proximity gap
27   inside methane.
28   Fixed cut on neutral particles in the StepManager.
29
30   Revision 1.6  2000/04/28 11:51:58  morsch
31    Dimensions of arrays hits and Ckov_data corrected.
32
33   Revision 1.5  2000/04/19 13:28:46  morsch
34   Major changes in geometry (parametrised), materials (updated) and
35   step manager (diagnostics) (JB, AM)
36
37 */
38
39
40
41 //////////////////////////////////////////////////////////
42 //  Manager and hits classes for set: RICH full version //
43 //////////////////////////////////////////////////////////
44
45 #include <TTUBE.h>
46 #include <TNode.h> 
47 #include <TRandom.h> 
48
49 #include "AliRICHv1.h"
50 #include "AliRun.h"
51 #include "AliMC.h"
52 #include "iostream.h"
53 #include "AliCallf77.h"
54 #include "AliConst.h" 
55 #include "AliPDG.h" 
56 #include "TGeant3.h"
57
58 ClassImp(AliRICHv1)
59     
60 //___________________________________________
61 AliRICHv1::AliRICHv1() : AliRICHv0()
62 {
63     //fChambers = 0;
64 }
65
66 //___________________________________________
67 AliRICHv1::AliRICHv1(const char *name, const char *title)
68     : AliRICHv0(name,title)
69 {
70     fCkov_number=0;
71     fFreon_prod=0;
72
73     fChambers = new TObjArray(7);
74     for (Int_t i=0; i<7; i++) {
75     
76         (*fChambers)[i] = new AliRICHChamber();  
77         
78     }
79 }
80
81
82 //___________________________________________
83 void AliRICHv1::StepManager()
84 {
85     Int_t          copy, id;
86     static Int_t   idvol;
87     static Int_t   vol[2];
88     Int_t          ipart;
89     static Float_t hits[18];
90     static Float_t Ckov_data[19];
91     TLorentzVector Position;
92     TLorentzVector Momentum;
93     Float_t        pos[3];
94     Float_t        mom[4];
95     Float_t        Localpos[3];
96     Float_t        Localmom[4];
97     Float_t        Localtheta,Localphi;
98     Float_t        theta,phi;
99     Float_t        destep, step;
100     Float_t        ranf[2];
101     Int_t          NPads;
102     Float_t        coscerenkov;
103     static Float_t eloss, xhit, yhit, tlength;
104     const  Float_t big=1.e10;
105        
106     TClonesArray &lhits = *fHits;
107     TGeant3 *geant3 = (TGeant3*) gMC;
108     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
109
110  //if (current->Energy()>1)
111    //{
112         
113     // Only gas gap inside chamber
114     // Tag chambers and record hits when track enters 
115     
116     idvol=-1;
117     id=gMC->CurrentVolID(copy);
118     Float_t cherenkov_loss=0;
119     //gAlice->KeepTrack(gAlice->CurrentTrack());
120     
121     gMC->TrackPosition(Position);
122     pos[0]=Position(0);
123     pos[1]=Position(1);
124     pos[2]=Position(2);
125     Ckov_data[1] = pos[0];                 // X-position for hit
126     Ckov_data[2] = pos[1];                 // Y-position for hit
127     Ckov_data[3] = pos[2];                 // Z-position for hit
128     //Ckov_data[11] = gAlice->CurrentTrack();
129
130     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH"); 
131     
132     /********************Store production parameters for Cerenkov photons************************/ 
133 //is it a Cerenkov photon? 
134     if (gMC->TrackPid() == 50000050) {          
135
136       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
137         //{                    
138           Float_t Ckov_energy = current->Energy();
139           //energy interval for tracking
140           if  (Ckov_energy > 5.6e-09 && Ckov_energy < 7.8e-09 )       
141             //if (Ckov_energy > 0)
142             {
143               if (gMC->IsTrackEntering()){                                     //is track entering?
144                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
145                   {                                                          //is it in freo?
146                     if (geant3->Gctrak()->nstep<1){                          //is it the first step?
147                       Int_t mother = current->GetFirstMother(); 
148                       
149                       //printf("Second Mother:%d\n",current->GetSecondMother());
150                       
151                       Ckov_data[10] = mother;
152                       Ckov_data[11] = gAlice->CurrentTrack();
153                       Ckov_data[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
154                       fCkov_number++;
155                       fFreon_prod=1;
156                       //printf("Index: %d\n",fCkov_number);
157                     }    //first step question
158                   }        //freo question
159                 
160                 if (geant3->Gctrak()->nstep<1){                                  //is it first step?
161                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
162                     {
163                       Ckov_data[12] = 2;
164                     }    //quarz question
165                 }        //first step question
166                 
167                 //printf("Before %d\n",fFreon_prod);
168               }   //track entering question
169               
170               if (Ckov_data[12] == 1)                                        //was it produced in Freon?
171                 //if (fFreon_prod == 1)
172                 {
173                   if (gMC->IsTrackEntering()){                                     //is track entering?
174                     //printf("Got in");
175                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?      
176                       {
177                         //printf("Got in\n");
178                         gMC->TrackMomentum(Momentum);
179                         mom[0]=Momentum(0);
180                         mom[1]=Momentum(1);
181                         mom[2]=Momentum(2);
182                         mom[3]=Momentum(3);
183                         // Z-position for hit
184                         
185                         
186                         /**************** Photons lost in second grid have to be calculated by hand************/ 
187                         
188                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
189                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
190                         gMC->Rndm(ranf, 1);
191                         //printf("grid calculation:%f\n",t);
192                         if (ranf[0] > t) {
193                           //geant3->StopTrack();
194                           Ckov_data[13] = 5;
195                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
196                           //printf("Lost one in grid\n");
197                         }
198                         /**********************************************************************************/
199                       }    //gap
200                     
201                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?      
202                       {
203                         gMC->TrackMomentum(Momentum);
204                         mom[0]=Momentum(0);
205                         mom[1]=Momentum(1);
206                         mom[2]=Momentum(2);
207                         mom[3]=Momentum(3);
208                         
209                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/ 
210                         /***********************Cerenkov phtons (always polarised)*************************/
211                         
212                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
213                         Float_t t = Fresnel(Ckov_energy*1e9,cophi,1);
214                         gMC->Rndm(ranf, 1);
215                         if (ranf[0] < t) {
216                           //geant3->StopTrack();
217                           Ckov_data[13] = 6;
218                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
219                           //printf("Lost by Fresnel\n");
220                         }
221                         /**********************************************************************************/
222                       }
223                   } //track entering?
224                   
225                   
226                   /********************Evaluation of losses************************/
227                   /******************still in the old fashion**********************/
228                   
229                   Int_t i1 = geant3->Gctrak()->nmec;            //number of physics mechanisms acting on the particle
230                   for (Int_t i = 0; i < i1; ++i) {
231                     //        Reflection loss 
232                     if (geant3->Gctrak()->lmec[i] == 106) {        //was it reflected
233                       Ckov_data[13]=10;
234                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
235                         Ckov_data[13]=1;
236                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
237                         Ckov_data[13]=2;
238                       //geant3->StopTrack();
239                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
240                     } //reflection question
241                     
242                     
243                     //        Absorption loss 
244                     else if (geant3->Gctrak()->lmec[i] == 101) {              //was it absorbed?
245                       Ckov_data[13]=20;
246                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
247                         Ckov_data[13]=11;
248                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
249                         Ckov_data[13]=12;
250                       if (gMC->CurrentVolID(copy) == gMC->VolId("META")) 
251                         Ckov_data[13]=13;
252                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP ")) 
253                         Ckov_data[13]=13;
254                       
255                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC")) 
256                         Ckov_data[13]=15;
257                       
258                       //        CsI inefficiency 
259                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
260                         Ckov_data[13]=16;
261                       }
262                       //geant3->StopTrack();
263                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
264                       //printf("Added cerenkov %d\n",fCkov_number);
265                     } //absorption question 
266                     
267                     
268                     //        Photon goes out of tracking scope 
269                     else if (geant3->Gctrak()->lmec[i] == 30) {                 //is it below energy treshold?
270                       Ckov_data[13]=21;
271                       //geant3->StopTrack();
272                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
273                     }   // energy treshold question         
274                   }  //number of mechanisms cycle
275                   /**********************End of evaluation************************/
276                 } //freon production question
277             } //energy interval question
278         //}//inside the proximity gap question
279     } //cerenkov photon question
280       
281     /**************************************End of Production Parameters Storing*********************/ 
282     
283     
284     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/ 
285     
286     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
287       //printf("Cerenkov\n");
288         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
289         {
290             
291           if (gMC->Edep() > 0.){
292                 gMC->TrackPosition(Position);
293                 gMC->TrackMomentum(Momentum);
294                 pos[0]=Position(0);
295                 pos[1]=Position(1);
296                 pos[2]=Position(2);
297                 mom[0]=Momentum(0);
298                 mom[1]=Momentum(1);
299                 mom[2]=Momentum(2);
300                 mom[3]=Momentum(3);
301                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
302                 Double_t rt = TMath::Sqrt(tc);
303                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
304                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
305                 gMC->Gmtod(pos,Localpos,1);                                                                    
306                 gMC->Gmtod(mom,Localmom,2);
307                 
308                 gMC->CurrentVolOffID(2,copy);
309                 vol[0]=copy;
310                 idvol=vol[0]-1;
311
312                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
313                         //->Sector(Localpos[0], Localpos[2]);
314                 //printf("Sector:%d\n",sector);
315
316                 /*if (gMC->TrackPid() == 50000051){
317                   fFeedbacks++;
318                   printf("Feedbacks:%d\n",fFeedbacks);
319                 }*/     
320                 
321                 ((AliRICHChamber*) (*fChambers)[idvol])
322                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
323                 if(idvol<7) {   
324                     Ckov_data[0] = gMC->TrackPid();        // particle type
325                     Ckov_data[1] = pos[0];                 // X-position for hit
326                     Ckov_data[2] = pos[1];                 // Y-position for hit
327                     Ckov_data[3] = pos[2];                 // Z-position for hit
328                     Ckov_data[4] = theta;                      // theta angle of incidence
329                     Ckov_data[5] = phi;                      // phi angle of incidence 
330                     Ckov_data[8] = (Float_t) fNPadHits;      // first padhit
331                     Ckov_data[9] = -1;                       // last pad hit
332                     Ckov_data[13] = 4;                       // photon was detected
333                     Ckov_data[14] = mom[0];
334                     Ckov_data[15] = mom[1];
335                     Ckov_data[16] = mom[2];
336                     
337                     destep = gMC->Edep();
338                     gMC->SetMaxStep(big);
339                     cherenkov_loss  += destep;
340                     Ckov_data[7]=cherenkov_loss;
341                     
342                     NPads = MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
343                     if (fNPadHits > (Int_t)Ckov_data[8]) {
344                         Ckov_data[8]= Ckov_data[8]+1;
345                         Ckov_data[9]= (Float_t) fNPadHits;
346                     }
347
348                     Ckov_data[17] = NPads;
349                     //printf("Npads:%d",NPads);
350                     
351                     //TClonesArray *Hits = RICH->Hits();
352                     AliRICHHit *mipHit =  (AliRICHHit*) (fHits->UncheckedAt(0));
353                     if (mipHit)
354                       {
355                         mom[0] = current->Px();
356                         mom[1] = current->Py();
357                         mom[2] = current->Pz();
358                         Float_t Mip_px = mipHit->fMomX;
359                         Float_t Mip_py = mipHit->fMomY;
360                         Float_t Mip_pz = mipHit->fMomZ;
361                         
362                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
363                         Float_t rt = TMath::Sqrt(r);
364                         Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;  
365                         Float_t Mip_rt = TMath::Sqrt(Mip_r);
366                         if ((rt*Mip_rt) > 0)
367                           {
368                             coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt);
369                           }
370                         else
371                           {
372                             coscerenkov = 0;
373                           }
374                         Float_t cherenkov = TMath::ACos(coscerenkov);
375                         Ckov_data[18]=cherenkov;
376                       }
377                     //if (sector != -1)
378                     //{
379                     AddHit(gAlice->CurrentTrack(),vol,Ckov_data);
380                     AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
381                     //}
382                 }
383             }
384         }
385     }
386     
387     /***********************************************End of photon hits*********************************************/
388     
389
390     /**********************************************Charged particles treatment*************************************/
391
392     else if (gMC->TrackCharge())
393     //else if (1 == 1)
394       {
395 //If MIP
396         /*if (gMC->IsTrackEntering())
397           {                
398             hits[13]=20;//is track entering?
399           }*/
400         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
401           {
402             fFreon_prod=1;
403           }
404
405         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
406 // Get current particle id (ipart), track position (pos)  and momentum (mom)
407             
408             gMC->CurrentVolOffID(3,copy);
409             vol[0]=copy;
410             idvol=vol[0]-1;
411
412             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
413                         //->Sector(Localpos[0], Localpos[2]);
414             //printf("Sector:%d\n",sector);
415             
416             gMC->TrackPosition(Position);
417             gMC->TrackMomentum(Momentum);
418             pos[0]=Position(0);
419             pos[1]=Position(1);
420             pos[2]=Position(2);
421             mom[0]=Momentum(0);
422             mom[1]=Momentum(1);
423             mom[2]=Momentum(2);
424             mom[3]=Momentum(3);
425             gMC->Gmtod(pos,Localpos,1);                                                                    
426             gMC->Gmtod(mom,Localmom,2);
427             
428             ipart  = gMC->TrackPid();
429             //
430             // momentum loss and steplength in last step
431             destep = gMC->Edep();
432             step   = gMC->TrackStep();
433   
434             //
435             // record hits when track enters ...
436             if( gMC->IsTrackEntering()) {
437 //              gMC->SetMaxStep(fMaxStepGas);
438                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
439                 Double_t rt = TMath::Sqrt(tc);
440                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
441                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
442                 
443
444                 Double_t Localtc = Localmom[0]*Localmom[0]+Localmom[2]*Localmom[2];
445                 Double_t Localrt = TMath::Sqrt(Localtc);
446                 Localtheta   = Float_t(TMath::ATan2(Localrt,Double_t(Localmom[1])))*kRaddeg;                       
447                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[2]),Double_t(Localmom[0])))*kRaddeg;    
448                 
449                 hits[0] = Float_t(ipart);         // particle type
450                 hits[1] = Localpos[0];                 // X-position for hit
451                 hits[2] = Localpos[1];                 // Y-position for hit
452                 hits[3] = Localpos[2];                 // Z-position for hit
453                 hits[4] = Localtheta;                  // theta angle of incidence
454                 hits[5] = Localphi;                    // phi angle of incidence 
455                 hits[8] = (Float_t) fNPadHits;    // first padhit
456                 hits[9] = -1;                     // last pad hit
457                 hits[13] = fFreon_prod;           // did id hit the freon?
458                 hits[14] = mom[0];
459                 hits[15] = mom[1];
460                 hits[16] = mom[2];
461
462                 tlength = 0;
463                 eloss   = 0;
464                 fFreon_prod = 0;
465         
466                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
467            
468                 
469                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
470                 xhit    = Localpos[0];
471                 yhit    = Localpos[2];
472                 // Only if not trigger chamber
473                 if(idvol<7) {
474                     //
475                     //  Initialize hit position (cursor) in the segmentation model 
476                     ((AliRICHChamber*) (*fChambers)[idvol])
477                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
478                 }
479             }
480             
481             // 
482             // Calculate the charge induced on a pad (disintegration) in case 
483             //
484             // Mip left chamber ...
485             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
486                 gMC->SetMaxStep(big);
487                 eloss   += destep;
488                 tlength += step;
489                 
490                                 
491                 // Only if not trigger chamber
492                 if(idvol<7) {
493                   if (eloss > 0) 
494                     {
495                       if(gMC->TrackPid() == kNeutron)
496                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
497                       NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
498                       hits[17] = NPads;
499                       //printf("Npads:%d",NPads);
500                     }
501                 }
502                 
503                 hits[6]=tlength;
504                 hits[7]=eloss;
505                 if (fNPadHits > (Int_t)hits[8]) {
506                     hits[8]= hits[8]+1;
507                     hits[9]= (Float_t) fNPadHits;
508                 }
509                 
510                 //if(sector !=-1)
511                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
512                 eloss = 0; 
513                 //
514                 // Check additional signal generation conditions 
515                 // defined by the segmentation
516                 // model (boundary crossing conditions) 
517             } else if 
518                 (((AliRICHChamber*) (*fChambers)[idvol])
519                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
520             {
521                 ((AliRICHChamber*) (*fChambers)[idvol])
522                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
523                 if (eloss > 0) 
524                   {
525                     if(gMC->TrackPid() == kNeutron)
526                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
527                     NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
528                     hits[17] = NPads;
529                     //printf("Npads:%d",NPads);
530                   }
531                 xhit     = Localpos[0];
532                 yhit     = Localpos[2]; 
533                 eloss    = destep;
534                 tlength += step ;
535                 //
536                 // nothing special  happened, add up energy loss
537             } else {        
538                 eloss   += destep;
539                 tlength += step ;
540             }
541         }
542       }
543     /*************************************************End of MIP treatment**************************************/
544    //}
545 }