ab82979e5c91d9a991f109e99299a291be88ea29
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterInput.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 #include <TClonesArray.h>
19 #include <TMinuit.h>
20
21 #include "AliRun.h"
22 #include "AliMUON.h"
23 #include "AliMUONChamber.h"
24 #include "AliMUONConstants.h"
25 #include "AliMUONClusterInput.h"
26 #include "AliMUONMathieson.h"
27 #include "AliMUONRawCluster.h"
28 #include "AliMUONDigit.h"
29 #include "AliLog.h"
30
31 ClassImp(AliMUONClusterInput)
32
33 AliMUONClusterInput* AliMUONClusterInput::fgClusterInput = 0; 
34 TMinuit* AliMUONClusterInput::fgMinuit = 0; 
35 AliMUONMathieson* AliMUONClusterInput::fgMathieson = 0; 
36
37 AliMUONClusterInput::AliMUONClusterInput()
38   : TObject(),
39     fCluster(0),
40     fChargeCorrel(1.),
41     fSegmentationType(1),
42     fDetElemId(0)
43   
44 {
45   fDigits[0]=0;
46   fDigits[1]=0;
47   fSegmentation[0]=0;
48   fSegmentation[1]=0;
49   fSegmentation2[0]=0;
50   fSegmentation2[1]=0;
51 }
52
53 AliMUONClusterInput* AliMUONClusterInput::Instance()
54 {
55 // return pointer to the singleton instance
56     if (fgClusterInput == 0) {
57         fgClusterInput = new AliMUONClusterInput();
58         fgMinuit = new TMinuit(8);
59     }
60     
61     return fgClusterInput;
62 }
63
64 AliMUONClusterInput::~AliMUONClusterInput()
65 {
66 // Destructor
67     delete fgMinuit;
68     delete fgMathieson;
69 }
70
71 AliMUONClusterInput::AliMUONClusterInput(const AliMUONClusterInput& clusterInput):TObject(clusterInput)
72 {
73 // Protected copy constructor
74
75   AliFatal("Not implemented.");
76 }
77
78 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig1, TClonesArray* dig2)
79 {
80 // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
81     fChamber=chamber;
82     fDigits[0]=dig1;
83     fDigits[1]=dig2; 
84     fNDigits[0]=dig1->GetEntriesFast();
85     fNDigits[1]=dig2->GetEntriesFast();
86     
87     AliMUON *pMUON;
88     AliMUONChamber* iChamber;
89
90     pMUON = (AliMUON*) gAlice->GetModule("MUON");
91     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
92       AliFatal("Wrong segmentation type");
93
94     iChamber =  &(pMUON->Chamber(chamber));
95     fgMathieson = new AliMUONMathieson();
96
97     fSegmentation[0]=iChamber->SegmentationModel(1);
98     fSegmentation[1]=iChamber->SegmentationModel(2);
99
100     fNseg = 2;
101     if (chamber < AliMUONConstants::NTrackingCh()) {
102       if (chamber > 1 ) {
103         fgMathieson->SetPitch(AliMUONConstants::Pitch());
104         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
105         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
106         fChargeCorrel = AliMUONConstants::ChargeCorrel();
107       } else {
108         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
109         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
110         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
111         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
112       }
113     }
114 }
115
116 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig)
117 {
118 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
119     fDigits[0]=dig;
120     AliMUON *pMUON;
121     AliMUONChamber* iChamber;
122
123     pMUON = (AliMUON*) gAlice->GetModule("MUON");
124     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
125       AliFatal("Wrong segmentation type");
126
127     iChamber =  &(pMUON->Chamber(chamber));
128     fSegmentation[0]=iChamber->SegmentationModel(1);
129
130     fNseg=1;
131 }
132 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig1, TClonesArray* dig2)
133 {
134   // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
135     fChamber = chamber;
136     fDetElemId = idDE;
137     fDigits[0]  = dig1;
138     fDigits[1]  = dig2; 
139     fNDigits[0] = dig1->GetEntriesFast();
140     fNDigits[1] = dig2->GetEntriesFast();
141     
142     AliMUON *pMUON;
143     AliMUONChamber* iChamber;
144
145     pMUON = (AliMUON*) gAlice->GetModule("MUON");
146     iChamber =  &(pMUON->Chamber(chamber));
147
148     fgMathieson = new AliMUONMathieson();
149     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
150       AliFatal("Wrong segmentation type");
151
152     fSegmentation2[0]=iChamber->SegmentationModel2(1);
153     fSegmentation2[1]=iChamber->SegmentationModel2(2);
154
155     fNseg = 2;
156     if (chamber < AliMUONConstants::NTrackingCh()) {
157       if (chamber > 1 ) {
158         fgMathieson->SetPitch(AliMUONConstants::Pitch());
159         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
160         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
161         fChargeCorrel = AliMUONConstants::ChargeCorrel();
162       } else {
163         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
164         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
165         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
166         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
167       }
168     }
169 }
170
171 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig)
172 {
173 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
174
175     fChamber = chamber;
176     fDetElemId = idDE;
177     fDigits[0] = dig;
178
179     AliMUON *pMUON;
180     AliMUONChamber* iChamber;
181
182     pMUON = (AliMUON*) gAlice->GetModule("MUON");
183     iChamber =  &(pMUON->Chamber(chamber));
184     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
185       AliFatal("Wrong segmentation type");
186     
187     fSegmentation2[0]=iChamber->SegmentationModel2(1);
188
189     fNseg=1;
190 }
191
192 void  AliMUONClusterInput::SetCluster(AliMUONRawCluster* cluster)
193 {
194 // Set the current cluster
195   //PH printf("\n %p \n", cluster);
196   fCluster=cluster;
197   Float_t qtot;
198   Int_t   i, cath, ix, iy;
199   AliMUONDigit* digit;
200   fNmul[0]=cluster->GetMultiplicity(0);
201   fNmul[1]=cluster->GetMultiplicity(1);
202   //PH printf("\n %p %p ", fDigits[0], fDigits[1]);
203   
204   for (cath=0; cath<2; cath++) {
205     qtot=0;
206     for (i=0; i<fNmul[cath]; i++) {
207       // pointer to digit
208       digit =(AliMUONDigit*)
209                 (fDigits[cath]->UncheckedAt(cluster->GetIndex(i,cath)));
210             // pad coordinates
211             ix = digit->PadX();
212             iy = digit->PadY();
213             // pad charge
214             fCharge[i][cath] = digit->Signal();
215             // pad centre coordinates
216 //          fSegmentation[cath]->GetPadCxy(ix, iy, x, y);
217             // globals kUsed in fitting functions
218             fix[i][cath]=ix;
219             fiy[i][cath]=iy;
220             // total charge per cluster
221             qtot+=fCharge[i][cath];
222             // Current z
223             Float_t xc, yc;
224             if (fSegmentationType == 1)
225               fSegmentation[cath]->GetPadC(ix,iy,xc,yc,fZ);
226             else 
227               fSegmentation2[cath]->GetPadC(fDetElemId,ix,iy,xc,yc,fZ);
228         } // loop over cluster digits
229         fQtot[cath]=qtot;
230         fChargeTot[cath]=Int_t(qtot);  
231     }  // loop over cathodes
232 }
233
234
235
236 Float_t AliMUONClusterInput::DiscrChargeS1(Int_t i,Double_t *par) 
237 {
238 // Compute the charge on first cathod only.
239 return DiscrChargeCombiS1(i,par,0);
240 }
241
242 Float_t AliMUONClusterInput::DiscrChargeCombiS1(Int_t i,Double_t *par, Int_t cath) 
243 {
244 // par[0]    x-position of cluster
245 // par[1]    y-position of cluster
246
247     Float_t q1;
248     if (fSegmentationType == 1) {
249
250       fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
251       //  First Cluster
252       fSegmentation[cath]->SetHit(par[0],par[1],fZ);
253       q1 = fgMathieson->IntXY(fSegmentation[cath]);
254
255     } else {
256
257       fSegmentation2[cath]->SetPad(fDetElemId, fix[i][cath], fiy[i][cath]);
258       //  First Cluster
259       fSegmentation2[cath]->SetHit(fDetElemId, par[0],par[1],fZ);
260       q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
261    }
262     
263    Float_t value = fQtot[cath]*q1;
264    return value;
265 }
266
267
268 Float_t AliMUONClusterInput::DiscrChargeS2(Int_t i,Double_t *par) 
269 {
270 // par[0]    x-position of first  cluster
271 // par[1]    y-position of first  cluster
272 // par[2]    x-position of second cluster
273 // par[3]    y-position of second cluster
274 // par[4]    charge fraction of first  cluster
275 // 1-par[4]  charge fraction of second cluster
276
277   Float_t q1, q2;
278   
279   if (fSegmentationType == 1) {
280
281     fSegmentation[0]->SetPad(fix[i][0], fiy[i][0]);
282     //  First Cluster
283     fSegmentation[0]->SetHit(par[0],par[1],fZ);
284     q1 = fgMathieson->IntXY(fSegmentation[0]);
285     
286     //  Second Cluster
287     fSegmentation[0]->SetHit(par[2],par[3],fZ);
288     q2 = fgMathieson->IntXY(fSegmentation[0]);
289
290   } else {
291
292     fSegmentation2[0]->SetPad(fDetElemId, fix[i][0], fiy[i][0]);
293     //  First Cluster
294     fSegmentation2[0]->SetHit(fDetElemId, par[0],par[1],fZ);
295     q1 = fgMathieson->IntXY(fSegmentation[0]);
296     
297     //  Second Cluster
298     fSegmentation2[0]->SetHit(fDetElemId,par[2],par[3],fZ);
299     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[0]);
300   }
301
302   Float_t value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
303   return value;
304 }
305
306 Float_t AliMUONClusterInput::DiscrChargeCombiS2(Int_t i,Double_t *par, Int_t cath) 
307 {
308 // par[0]    x-position of first  cluster
309 // par[1]    y-position of first  cluster
310 // par[2]    x-position of second cluster
311 // par[3]    y-position of second cluster
312 // par[4]    charge fraction of first  cluster - first cathode
313 // 1-par[4]  charge fraction of second cluster 
314 // par[5]    charge fraction of first  cluster - second cathode
315
316   Float_t q1, q2;
317
318   if (fSegmentationType == 1) {
319
320     fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
321     //  First Cluster
322     fSegmentation[cath]->SetHit(par[0],par[1],fZ);
323     q1 = fgMathieson->IntXY(fSegmentation[cath]);
324     
325     //  Second Cluster
326     fSegmentation[cath]->SetHit(par[2],par[3],fZ);
327     q2 = fgMathieson->IntXY(fSegmentation[cath]);
328
329   } else {
330     fSegmentation2[cath]->SetPad(fDetElemId,fix[i][cath], fiy[i][cath]);
331     //  First Cluster
332     fSegmentation2[cath]->SetHit(fDetElemId,par[0],par[1],fZ);
333     q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
334     
335     //  Second Cluster
336     fSegmentation2[cath]->SetHit(fDetElemId,par[2],par[3],fZ);
337     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
338   }
339   Float_t value;
340   if (cath==0) {
341     value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
342   } else {
343     value = fQtot[1]*(par[5]*q1+(1.-par[5])*q2);
344   }
345   return value;
346 }
347
348 AliMUONClusterInput& AliMUONClusterInput
349 ::operator = (const AliMUONClusterInput& rhs)
350 {
351 // Protected assignement operator
352
353   if (this == &rhs) return *this;
354
355   AliFatal("Not implemented.");
356     
357   return *this;  
358 }