b868de8241b359e92c5e89f3050a8b4007ca67b8
[u/mrichter/AliRoot.git] / MUON / AliMUONChamber.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 // --- ROOT includes ---
19 #include <TRandom.h>
20 #include <TMath.h>
21 #include "AliRun.h"
22
23
24 // --- MUON includes ---
25 #include "AliMUON.h"
26 #include "AliMUONChamber.h"
27 #include "AliMUONGeometryModule.h"
28 #include "AliMUONHit.h"
29 #include "AliLog.h"
30
31 ClassImp(AliMUONChamber)        
32
33 AliMUONChamber::AliMUONChamber()
34   : TObject(), 
35     fId(0),
36     fdGas(0.),
37     fdAlu(0.),
38     fZ(0.),
39     fnsec(1),
40     frMin(0.),
41     frMax(0.),
42     fCurrentCorrel(1), // to avoid mistakes if ChargeCorrelInit is not called
43     fSegmentation(0),
44     fSegmentation2(0),
45     fResponse(0),
46     fGeometry(0),
47     fMUON(0)
48 {
49 // Default constructor
50 }
51
52 //_______________________________________________________
53 AliMUONChamber::AliMUONChamber(Int_t id) 
54   : TObject(), 
55     fId(id),
56     fdGas(0.),
57     fdAlu(0.),
58     fZ(0.),
59     fnsec(1),
60     frMin(0.),
61     frMax(0.),
62     fCurrentCorrel(1), // to avoid mistakes if ChargeCorrelInit is not called
63     fSegmentation(0),
64     fSegmentation2(0),
65     fResponse(0),
66     fGeometry(0),
67     fMUON(0)
68 {
69
70     // muon
71     fMUON = (AliMUON*)gAlice->GetModule("MUON");
72     if (!fMUON) {
73       AliFatal("MUON detector not defined.");
74       return;
75     }  
76 // Construtor with chamber id 
77     fSegmentation = new TObjArray(2);
78     fSegmentation->AddAt(0,0);
79     fSegmentation->AddAt(0,1);
80     
81     // new segmentation
82     fSegmentation2 = new TObjArray(2);
83     fSegmentation2->AddAt(0,0);
84     fSegmentation2->AddAt(0,1);
85  
86     fGeometry = new AliMUONGeometryModule(fId);
87
88 }
89
90 //_______________________________________________________
91 AliMUONChamber::AliMUONChamber(const AliMUONChamber& rChamber)
92   : TObject(rChamber)
93 {
94   // Protected copy constructor
95
96   AliFatal("Not implemented.");
97   // Dummy copy constructor
98 }
99
100 //_______________________________________________________
101 AliMUONChamber::~AliMUONChamber() 
102 {
103   // Destructor
104   if (fMUON->WhichSegmentation() == 1) {
105     if (fSegmentation) {
106       fSegmentation->Delete();
107       delete fSegmentation;
108     } 
109   } else {
110     if (fSegmentation2) {
111       fSegmentation2->Delete();
112       delete fSegmentation2;
113     }
114   }
115   delete fGeometry;
116 }
117
118 //_______________________________________________________
119 AliMUONChamber & AliMUONChamber::operator =(const AliMUONChamber& rhs)
120 {
121   // Protected assignement operator
122
123   if (this == &rhs) return *this;
124
125   AliFatal("Not implemented.");
126     
127   return *this;  
128 }
129
130 //_______________________________________________________
131 Bool_t  AliMUONChamber::IsSensId(Int_t volId) const 
132 {
133   // Returns true if the volume specified by volId is in the list
134   // of sesitive volumes for this chamber
135
136   return fGeometry->IsSensitiveVolume(volId);
137 }  
138
139 //_______________________________________________________
140 void AliMUONChamber::Init()
141 {
142   // Initalisation ..
143   //
144   // ... for chamber segmentation
145
146   if (fSegmentation->At(0)) 
147     ((AliSegmentation *) fSegmentation->At(0))->Init(fId);
148
149   if (fnsec==2) {
150     if (fSegmentation->At(1))
151       ((AliSegmentation *) fSegmentation->At(1))->Init(fId);
152   }
153
154   
155 }
156
157 //_________________________________________________________________
158 Int_t   AliMUONChamber::SigGenCond(Float_t x, Float_t y, Float_t z)
159 {
160   // Ask segmentation if signal should be generated 
161
162   if (fnsec==1) {
163     return ((AliSegmentation*) fSegmentation->At(0))
164       ->SigGenCond(x, y, z) ;
165   } else {
166     return (((AliSegmentation*) fSegmentation->At(0))
167             ->SigGenCond(x, y, z)) ||
168       (((AliSegmentation*) fSegmentation->At(1))
169        ->SigGenCond(x, y, z)) ;
170   }
171   
172 }
173
174 //_________________________________________________________________
175 void    AliMUONChamber::SigGenInit(Float_t x, Float_t y, Float_t z)
176 {
177   //
178   // Initialisation of segmentation for hit
179   //  
180
181
182   if (fnsec==1) {
183     ((AliSegmentation*) fSegmentation->At(0))->SigGenInit(x, y, z) ;
184   } else {
185     ((AliSegmentation*) fSegmentation->At(0))->SigGenInit(x, y, z) ;
186     ((AliSegmentation*) fSegmentation->At(1))->SigGenInit(x, y, z) ;
187   }
188   
189 }
190
191 //_____________________________________________________
192 void AliMUONChamber::ChargeCorrelationInit() {
193   // Initialisation of charge correlation for current hit
194   // the value is stored, and then used by Disintegration
195   if (fnsec==1) 
196     fCurrentCorrel =1;
197   else 
198     // exponential is here to avoid eventual problems in 0 
199     // factor 2 because chargecorrel is q1/q2 and not q1/qtrue
200     fCurrentCorrel = TMath::Exp(gRandom->Gaus(0,fResponse->ChargeCorrel()/2));
201 }
202
203 //_______________________________________________________
204 void AliMUONChamber::DisIntegration(Float_t eloss, Float_t /*tof*/, 
205                                     Float_t xhit, Float_t yhit, Float_t zhit,
206                                     Int_t& nnew,Float_t newclust[6][500]) 
207 {
208   //    
209   //  Generates pad hits (simulated cluster) 
210   //  using the segmentation and the response model 
211   Float_t dx, dy;
212   //
213   // Width of the integration area
214   //
215   dx=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
216   dy=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
217   //
218   // Get pulse height from energy loss
219   Float_t qtot = fResponse->IntPH(eloss);
220   //
221   // Loop Over Pads
222     
223   Float_t qp; 
224   nnew=0;
225     
226   // Cathode plane loop
227   for (Int_t i=1; i<=fnsec; i++) {
228     Float_t qcath = qtot * (i==1? fCurrentCorrel : 1/fCurrentCorrel);
229     AliSegmentation* segmentation=
230       (AliSegmentation*) fSegmentation->At(i-1);
231     for (segmentation->FirstPad(xhit, yhit, zhit, dx, dy); 
232          segmentation->MorePads(); 
233          segmentation->NextPad()) 
234       {
235         qp=fResponse->IntXY(segmentation);
236         qp=TMath::Abs(qp);
237         //
238         //
239         if (qp > 1.e-4) 
240           {
241             if (nnew >= 500) // Perform a bounds check on nnew since it is assumed
242               // newclust only contains 500 elements.
243               {
244                 AliError("Limit of 500 pad responses reached.");
245                 return;
246               };
247             //
248             // --- store signal information
249             newclust[0][nnew]=qcath;                     // total charge
250             newclust[1][nnew]=segmentation->Ix();       // ix-position of pad
251             newclust[2][nnew]=segmentation->Iy();       // iy-position of pad
252             newclust[3][nnew]=qp * qcath;                // charge on pad
253             newclust[4][nnew]=segmentation->ISector();  // sector id
254             newclust[5][nnew]=(Float_t) i;              // counter
255             nnew++;
256                 
257           }
258       } // Pad loop
259   } // Cathode plane loop
260 }
261
262 //_______________________________________________________
263 void AliMUONChamber::InitGeo(Float_t /*zpos*/)
264 {
265   //    sensitive gas gap
266   fdGas= 0.5;
267   //    3% radiation length of aluminum (X0=8.9 cm)      
268   fdAlu= 3.0/100*8.9;
269 }
270 //_______________________________________________________
271 //
272 //                  NEW SEGMENTATION
273 //_______________________________________________________
274 void AliMUONChamber::Init(Int_t flag)
275 {
276   // Initalisation ..
277   //
278   // ... for chamber segmentation
279
280   if (!flag)    AliFatal("wrong segmentation type.");
281
282
283   if (fSegmentation2->At(0)) 
284     ((AliMUONGeometrySegmentation*) fSegmentation2->At(0))->Init(fId);
285
286   if (fnsec==2) {
287     if (fSegmentation2->At(1))
288       ((AliMUONGeometrySegmentation*) fSegmentation2->At(1))->Init(fId);
289   }
290 }
291
292 //_________________________________________________________________
293 Int_t   AliMUONChamber::SigGenCond(AliMUONHit *hit)
294 {
295   // Ask segmentation if signal should be generated 
296
297   Float_t x = hit->X();
298   Float_t y = hit->Y();
299   Float_t z = hit->Z();
300   Int_t  id = hit->DetElemId();
301   
302   if (fnsec==1) {
303     return  ((AliMUONGeometrySegmentation*)fSegmentation2->At(0))->SigGenCond(id, x, y, z);
304   } else {
305     return (((AliMUONGeometrySegmentation*) fSegmentation2->At(0))
306             ->SigGenCond(id, x, y, z)) ||
307       (((AliMUONGeometrySegmentation*) fSegmentation2->At(1))
308        ->SigGenCond(id, x, y, z)) ;
309   }
310 }
311
312 //_________________________________________________________________
313 void    AliMUONChamber::SigGenInit(AliMUONHit *hit)
314 {
315   //
316   // Initialisation of segmentation for hit
317   //  
318   Float_t x = hit->X();
319   Float_t y = hit->Y();
320   Float_t z = hit->Z();
321   Int_t  id = hit->DetElemId();
322
323   if (fnsec==1) {
324     ((AliMUONGeometrySegmentation*) fSegmentation2->At(0))->SigGenInit(id, x, y, z) ;
325   } else {
326     ((AliMUONGeometrySegmentation*) fSegmentation2->At(0))->SigGenInit(id, x, y, z) ;
327     ((AliMUONGeometrySegmentation*) fSegmentation2->At(1))->SigGenInit(id, x, y, z) ;
328   }
329 }
330
331 //_______________________________________________________
332 void AliMUONChamber::DisIntegration(AliMUONHit *hit, 
333                                     Int_t& nnew,Float_t newclust[6][500]) 
334 {
335   //    
336   //  Generates pad hits (simulated cluster) 
337   //  using the segmentation and the response model 
338   Float_t dx, dy;
339
340   Float_t  xhit = hit->X();
341   Float_t  yhit = hit->Y();
342   Float_t  zhit = hit->Z();
343   Int_t      id = hit->DetElemId();
344   Float_t eloss = hit->Eloss();
345
346   //
347   // Width of the integration area
348   //
349   dx=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
350   dy=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
351   //
352   // Get pulse height from energy loss
353   Float_t qtot = fResponse->IntPH(eloss);
354   //
355   // Loop Over Pads
356     
357   Float_t qp; 
358   nnew=0;
359     
360   // Cathode plane loop
361   for (Int_t i=1; i<=fnsec; i++) {
362     Float_t qcath = qtot * (i==1? fCurrentCorrel : 1/fCurrentCorrel);
363
364     AliMUONGeometrySegmentation* segmentation=
365       (AliMUONGeometrySegmentation*) fSegmentation2->At(i-1);
366
367     for (segmentation->FirstPad(id, xhit, yhit, zhit, dx, dy); 
368          segmentation->MorePads(id); 
369          segmentation->NextPad(id)) 
370       {
371         qp=fResponse->IntXY(id, segmentation);
372         qp=TMath::Abs(qp);
373         //
374         //
375         if (qp > 1.e-4) 
376           {
377             if (nnew >= 500) // Perform a bounds check on nnew since it is assumed
378               // newclust only contains 500 elements.
379               {
380                 AliError("Limit of 500 pad responses reached.");
381                 return;
382               };
383             //
384             // --- store signal information
385             newclust[0][nnew]=qcath;                     // total charge
386             newclust[1][nnew]=segmentation->Ix();       // ix-position of pad
387             newclust[2][nnew]=segmentation->Iy();       // iy-position of pad
388             newclust[3][nnew]=qp * qcath;                // charge on pad
389             newclust[4][nnew]=segmentation->ISector();  // sector id
390             newclust[5][nnew]=(Float_t) i;              // counter
391             nnew++;
392                 
393           }
394       } // Pad loop
395   } // Cathode plane loop
396 }