d60b27812752bd37a3ff74ff3439fd756adbd7fc
[u/mrichter/AliRoot.git] / MUON / AliMUON.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 $Log$
17 Revision 1.25  2000/06/29 12:34:09  morsch
18 AliMUONSegmentation class has been made independent of AliMUONChamber. This makes
19 it usable with any other geometry class. The link to the object to which it belongs is
20 established via an index. This assumes that there exists a global geometry manager
21 from which the pointer to the parent object can be obtained (in our case gAlice).
22
23 Revision 1.24  2000/06/28 15:16:35  morsch
24 (1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
25 to allow development of slat-muon chamber simulation and reconstruction code in the MUON
26 framework. The changes should have no side effects (mostly dummy arguments).
27 (2) Hit disintegration uses 3-dim hit coordinates to allow simulation
28 of chambers with overlapping modules (MakePadHits, Disintegration).
29
30 Revision 1.23  2000/06/28 12:19:17  morsch
31 More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
32 cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
33 AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
34 It requires two cathode planes. Small modifications in the code will make it usable for
35 one cathode plane and, hence, more general (for test beam data).
36 AliMUONClusterFinder is now obsolete.
37
38 Revision 1.22  2000/06/28 08:06:10  morsch
39 Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
40 algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
41 It also naturally takes care of the TMinuit instance.
42
43 Revision 1.21  2000/06/27 08:54:41  morsch
44 Problems with on constant array sizes (in hitMap, nmuon, xhit, yhit) corrected.
45
46 Revision 1.20  2000/06/26 14:02:38  morsch
47 Add class AliMUONConstants with MUON specific constants using static memeber data and access methods.
48
49 Revision 1.19  2000/06/22 13:40:51  morsch
50 scope problem on HP, "i" declared once
51 pow changed to TMath::Power (PH, AM)
52
53 Revision 1.18  2000/06/15 07:58:48  morsch
54 Code from MUON-dev joined
55
56 Revision 1.14.4.17  2000/06/14 14:36:46  morsch
57 - add TriggerCircuit (PC)
58 - add GlobalTrigger and LocalTrigger and specific methods (PC)
59
60 Revision 1.14.4.16  2000/06/09 21:20:28  morsch
61 Most coding rule violations corrected
62
63 Revision 1.14.4.15  2000/05/02 09:54:32  morsch
64 RULE RN17 violations corrected
65
66 Revision 1.14.4.12  2000/04/26 12:25:02  morsch
67 Code revised by P. Crochet:
68 - Z position of TriggerChamber changed according to A.Tournaire Priv.Comm.
69 - ToF included in the method MakePadHits
70 - inner radius of flange between beam shielding and trigger corrected
71 - Trigger global volume updated (according to the new geometry)
72
73 Revision 1.14.4.11  2000/04/19 19:42:08  morsch
74 Some changes of variable names curing viols and methods concerning
75 correlated clusters removed.
76
77 Revision 1.14.4.10  2000/03/22 16:44:07  gosset
78 Memory leak suppressed in function Digitise:
79 p_adr->Delete() instead of Clear (I.Chevrot and A.Baldisseri)
80
81 Revision 1.14.4.9  2000/03/20 18:15:25  morsch
82 Positions of trigger chambers corrected (P.C.)
83
84 Revision 1.14.4.8  2000/02/21 15:38:01  morsch
85 Call to AddHitList introduced to make this version compatible with head.
86
87 Revision 1.14.4.7  2000/02/20 07:45:53  morsch
88 Bugs in Trigger part of BuildGeomemetry corrected (P.C)
89
90 Revision 1.14.4.6  2000/02/17 14:28:54  morsch
91 Trigger included into initialization and digitization
92
93 Revision 1.14.4.5  2000/02/15 10:02:58  morsch
94 Log messages of previous revisions added
95
96 Revision 1.14.4.2  2000/02/04 10:57:34  gosset
97 Z position of the chambers:
98 it was the Z position of the stations;
99 it is now really the Z position of the chambers.
100    !!!! WARNING: THE CALLS TO "AliMUONChamber::SetZPOS"
101    !!!!                   AND "AliMUONChamber::ZPosition"
102    !!!! HAVE TO BE CHANGED TO "AliMUONChamber::"SetZ"
103    !!!!                   AND "AliMUONChamber::Z"
104
105 Revision 1.14.4.3  2000/02/04 16:19:04  gosset
106 Correction for mis-spelling of NCH 
107
108 Revision 1.14.4.4  2000/02/15 09:43:38  morsch
109 Log message added
110
111 */
112
113
114 ///////////////////////////////////////////////
115 //  Manager and hits classes for set:MUON     //
116 ////////////////////////////////////////////////
117
118 #include <TTUBE.h>
119 #include <TBRIK.h>
120 #include <TRotMatrix.h>
121 #include <TNode.h> 
122 #include <TTree.h> 
123 #include <TRandom.h> 
124 #include <TObject.h>
125 #include <TVector.h>
126 #include <TObjArray.h>
127 #include <TMinuit.h>
128 #include <TParticle.h>
129 #include <TROOT.h>
130 #include <TFile.h>
131 #include <TNtuple.h>
132 #include <TCanvas.h>
133 #include <TPad.h>
134 #include <TDirectory.h>
135 #include <TObjectTable.h>
136 #include <AliPDG.h>
137 #include <TTUBE.h>
138
139 #include "AliMUON.h"
140 #include "AliMUONHit.h"
141 #include "AliMUONPadHit.h"
142 #include "AliMUONDigit.h"
143 #include "AliMUONTransientDigit.h"
144 #include "AliMUONRawCluster.h"
145 #include "AliMUONLocalTrigger.h"
146 #include "AliMUONGlobalTrigger.h"
147 #include "AliHitMap.h"
148 #include "AliMUONHitMapA1.h"
149 #include "AliMUONChamberTrigger.h"
150 #include "AliMUONConstants.h"
151 #include "AliMUONClusterFinderVS.h"
152 #include "AliMUONTriggerDecision.h"
153 #include "AliRun.h"
154 #include "AliMC.h"
155 #include "AliMUONClusterInput.h"
156 #include "iostream.h"
157 #include "AliCallf77.h" 
158 #include "AliConst.h" 
159
160 // Defaults parameters for Z positions of chambers
161 // taken from values for "stations" in AliMUON::AliMUON
162 //     const Float_t zch[7]={528, 690., 975., 1249., 1449., 1610, 1710.};
163 // and from array "dstation" in AliMUONv1::CreateGeometry
164 //          Float_t dstation[5]={20., 20., 20, 20., 20.};
165 //     for tracking chambers,
166 //          according to (Z1 = zch - dstation) and  (Z2 = zch + dstation)
167 //          for the first and second chambers in the station, respectively,
168 // and from "DTPLANES" in AliMUONv1::CreateGeometry
169 //           const Float_t DTPLANES = 15.;
170 //     for trigger chambers,
171 //          according to (Z1 = zch) and  (Z2 = zch + DTPLANES)
172 //          for the first and second chambers in the station, respectively
173
174 ClassImp(AliMUON)
175 //___________________________________________
176 AliMUON::AliMUON()
177 {
178    fIshunt       = 0;
179    fHits         = 0;
180    fPadHits      = 0;
181    fNPadHits     = 0;
182    fDchambers    = 0;
183    fTriggerCircuits = 0;     // cp new design of AliMUONTriggerDecision
184    fNdch         = 0;
185    fRawClusters  = 0;
186    fNrawch       = 0;
187    fGlobalTrigger   = 0;
188    fNLocalTrigger   = 0;
189    fLocalTrigger    = 0;
190    fNLocalTrigger   = 0;
191 }
192  
193 //___________________________________________
194 AliMUON::AliMUON(const char *name, const char *title)
195        : AliDetector(name,title)
196 {
197 //Begin_Html
198 /*
199 <img src="gif/alimuon.gif">
200 */
201 //End_Html
202
203    fHits     = new TClonesArray("AliMUONHit",1000);
204    gAlice->AddHitList(fHits);
205    fPadHits = new TClonesArray("AliMUONPadHit",10000);
206    fNPadHits  =  0;
207    fIshunt     =  0;
208
209    fNdch      = new Int_t[AliMUONConstants::NCh()];
210
211    fDchambers = new TObjArray(AliMUONConstants::NCh());
212
213    Int_t i;
214    
215    for (i=0; i<AliMUONConstants::NCh() ;i++) {
216        (*fDchambers)[i] = new TClonesArray("AliMUONDigit",10000); 
217        fNdch[i]=0;
218    }
219
220    fNrawch      = new Int_t[AliMUONConstants::NTrackingCh()];
221
222    fRawClusters = new TObjArray(AliMUONConstants::NTrackingCh());
223
224    for (i=0; i<AliMUONConstants::NTrackingCh();i++) {
225        (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); 
226        fNrawch[i]=0;
227    }
228    cout << " here " << "\n";
229
230    fGlobalTrigger = new TClonesArray("AliMUONGlobalTrigger",1);    
231    fNGlobalTrigger = 0;
232    fLocalTrigger  = new TClonesArray("AliMUONLocalTrigger",234);    
233    fNLocalTrigger = 0;
234
235 //   
236 // Transport angular cut
237    fAccCut=0;
238    fAccMin=2;
239    fAccMax=9;
240
241    SetMarkerColor(kRed);
242 //
243 //
244 //
245 //
246
247     Int_t ch;
248
249     fChambers = new TObjArray(AliMUONConstants::NCh());
250
251     // Loop over stations
252     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
253       // Loop over 2 chambers in the station
254         for (Int_t stCH = 0; stCH < 2; stCH++) {
255 //
256 //    
257 //    Default Parameters for Muon Tracking Stations
258
259
260             ch = 2 * st + stCH;
261 //
262             if (ch < AliMUONConstants::NTrackingCh()) {
263                 (*fChambers)[ch] = new AliMUONChamber(ch);
264             } else {
265                 (*fChambers)[ch] = new AliMUONChamberTrigger(ch);
266             }
267             
268             AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[ch];
269             
270             chamber->SetGid(0);
271             // Default values for Z of chambers
272             chamber->SetZ(AliMUONConstants::DefaultChamberZ(ch));
273 //
274             chamber->InitGeo(AliMUONConstants::DefaultChamberZ(ch));
275             chamber->SetRInner(AliMUONConstants::Dmin(st)/2);
276             chamber->SetROuter(AliMUONConstants::Dmax(st)/2);
277 //
278         } // Chamber stCH (0, 1) in 
279     }     // Station st (0...)
280     fMaxStepGas=0.01; 
281     fMaxStepAlu=0.1; 
282     fMaxDestepGas=-1;
283     fMaxDestepAlu=-1;
284 //
285    fMaxIterPad   = 0;
286    fCurIterPad   = 0;
287
288    // cp new design of AliMUONTriggerDecision
289    fTriggerCircuits = new TObjArray(AliMUONConstants::NTriggerCircuit());
290    for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
291      (*fTriggerCircuits)[circ] = new AliMUONTriggerCircuit();     
292    }
293    // cp new design of AliMUONTriggerDecision
294
295 }
296  
297 //___________________________________________
298 AliMUON::AliMUON(const AliMUON& rMUON)
299 {
300 // Dummy copy constructor
301     ;
302     
303 }
304
305 AliMUON::~AliMUON()
306 {
307     printf("Calling AliMUON destructor !!!\n");
308     
309     Int_t i;
310     fIshunt  = 0;
311     delete fHits;
312     delete fPadHits;
313     
314     delete fGlobalTrigger;
315     fNGlobalTrigger = 0;
316     
317     delete fLocalTrigger;
318     fNLocalTrigger = 0;
319
320     for (i=0;i<AliMUONConstants::NCh();i++) {
321         delete (*fDchambers)[i];
322         fNdch[i]=0;
323     }
324     delete fDchambers;
325     
326     for (i=0;i<AliMUONConstants::NTrackingCh();i++) {
327         delete (*fRawClusters)[i];
328         fNrawch[i]=0;
329     }
330     delete fRawClusters;
331     
332     for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
333         delete (*fTriggerCircuits)[circ];
334     }
335     delete fTriggerCircuits;
336 }
337  
338 //___________________________________________
339 void AliMUON::AddHit(Int_t track, Int_t *vol, Float_t *hits)
340 {
341   TClonesArray &lhits = *fHits;
342   new(lhits[fNhits++]) AliMUONHit(fIshunt,track,vol,hits);
343 }
344 //___________________________________________
345 void AliMUON::AddPadHit(Int_t *clhits)
346 {
347    TClonesArray &lclusters = *fPadHits;
348    new(lclusters[fNPadHits++]) AliMUONPadHit(clhits);
349 }
350 //_____________________________________________________________________________
351 void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
352 {
353     //
354     // Add a MUON digit to the list
355     //
356
357     TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
358     new(ldigits[fNdch[id]++]) AliMUONDigit(tracks,charges,digits);
359 }
360
361 //_____________________________________________________________________________
362 void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
363 {
364     //
365     // Add a MUON digit to the list
366     //
367
368     TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
369     new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
370 }
371
372 //___________________________________________
373 void AliMUON::AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus,
374                                Int_t *singleUndef,
375                                Int_t *pairUnlike, Int_t *pairLike)
376 {
377 // add a MUON Global Trigger to the list (only one GlobalTrigger per event !)
378   TClonesArray &globalTrigger = *fGlobalTrigger;
379   new(globalTrigger[fNGlobalTrigger++]) 
380     AliMUONGlobalTrigger(singlePlus, singleMinus,  singleUndef, pairUnlike, 
381                          pairLike);
382 }
383 //___________________________________________
384 void AliMUON::AddLocalTrigger(Int_t *localtr)
385 {
386 // add a MUON Local Trigger to the list
387   TClonesArray &localTrigger = *fLocalTrigger;
388   new(localTrigger[fNLocalTrigger++]) AliMUONLocalTrigger(localtr);
389 }
390
391 //___________________________________________
392 void AliMUON::BuildGeometry()
393 {
394     TNode *node, *nodeF, *top, *nodeS;
395     const int kColorMUON  = kBlue;
396     const int kColorMUON2 = kGreen; 
397     //
398     top=gAlice->GetGeometry()->GetNode("alice");
399 // MUON
400 //
401 //  z-Positions of Chambers
402     const Float_t kCz[7]={511., 686., 971., 1245., 1445., 1600, 1700.};
403 //  inner diameter (Xlenght for trigger chamber -> active area)
404     const Float_t kDmin[7]={ 35.,  47.,  67.,   86.,  100., 544., 544.};
405 //  outer diameter (Ylenght for trigger chamber -> active area)
406     const Float_t kDmax[7]={183., 245., 346.,  520.,  520., 612., 612.};
407
408     TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90,  0, 90, 90, 0, 0);
409     TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
410     TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
411     TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90,  0, 0, 0);
412     
413     Float_t rmin, rmax, dx, dy, dz, dr, xpos, ypos, zpos;
414     Float_t dzc1=4.;           // tracking chambers
415     Float_t dzc2=15.;          // trigger chambers
416     Float_t hole=102.;          // x-y hole around beam pipe for trig. chambers
417     Float_t zscale;            // scaling parameter trigger chambers
418     Float_t halfx, halfy;   
419     char nameChamber[9], nameSense[9], nameFrame[9], nameNode[8];
420     char nameSense1[9], nameSense2[9];    
421     for (Int_t i=0; i<7; i++) {
422         for (Int_t j=0; j<2; j++) {
423             Int_t id=2*i+j+1;
424             if (i<5) {               // tracking chambers
425                 if (j==0) {
426                     zpos=kCz[i]-dzc1;
427                 } else {
428                     zpos=kCz[i]+dzc1;
429                 }
430             } else {
431                 if (j==0) {
432                     zpos=kCz[i];
433                 } else {            
434                     zpos=kCz[i]+dzc2;
435                 }
436             }
437             sprintf(nameChamber,"C_MUON%d",id);
438             sprintf(nameSense,"S_MUON%d",id);
439             sprintf(nameSense1,"S1_MUON%d",id);
440             sprintf(nameSense2,"S2_MUON%d",id);
441             sprintf(nameFrame,"F_MUON%d",id);   
442             if (i<5) {                        // tracking chambers
443                 rmin = kDmin[i]/2.-3;
444                 rmax = kDmax[i]/2.+3;
445                 new TTUBE(nameChamber,"Mother","void",rmin,rmax,0.25,1.);
446                 rmin = kDmin[i]/2.;
447                 rmax = kDmax[i]/2.;
448                 new TTUBE(nameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
449                 dx=(rmax-rmin)/2;
450                 dy=3.;
451                 dz=0.25;
452                 TBRIK* frMUON = new TBRIK(nameFrame,"Frame","void",dx,dy,dz);
453                 top->cd();
454                 sprintf(nameNode,"MUON%d",100+id);
455                 node = new TNode(nameNode,"ChamberNode",nameChamber,0,0,zpos,"");
456                 node->SetLineColor(kColorMUON);
457                 fNodes->Add(node);
458                 node->cd();
459                 sprintf(nameNode,"MUON%d",200+id);
460                 node = new TNode(nameNode,"Sens. Region Node",nameSense,0,0,0,"");
461                 node->SetLineColor(kColorMUON);
462                 node->cd();
463                 dr=dx+rmin;
464                 sprintf(nameNode,"MUON%d",300+id);
465                 nodeF = new TNode(nameNode,"Frame0",frMUON,dr, 0, 0,rot000,"");
466                 nodeF->SetLineColor(kColorMUON);
467                 node->cd();
468                 sprintf(nameNode,"MUON%d",400+id);
469                 nodeF = new TNode(nameNode,"Frame1",frMUON,0 ,dr,0,rot090,"");
470                 nodeF->SetLineColor(kColorMUON);
471                 node->cd();
472                 sprintf(nameNode,"MUON%d",500+id);
473                 nodeF = new TNode(nameNode,"Frame2",frMUON,-dr,0,0,rot180,"");
474                 nodeF->SetLineColor(kColorMUON);
475                 node  ->cd();
476                 sprintf(nameNode,"MUON%d",600+id);   
477                 nodeF = new TNode(nameNode,"Frame3",frMUON,0,-dr,0,rot270,"");
478                 nodeF->SetLineColor(kColorMUON);
479             } else { 
480                 zscale=zpos/kCz[5];
481                 Float_t xsize=kDmin[i]*zscale;
482                 Float_t ysize=kDmax[i]*zscale;
483                 Float_t holeScaled=hole*zscale;
484                 
485                 halfx=xsize/2.+3.;
486                 halfy=ysize/2.+3.;          
487                 new TBRIK(nameChamber,"Mother","void",halfx,halfy,0.25);
488                 top->cd();
489                 sprintf(nameNode,"MUON%d",100+id);
490                 node = new TNode(nameNode,"Chambernode",nameChamber,0,0,zpos,"");
491                 node->SetLineColor(kColorMUON2);
492                 fNodes->Add(node);
493                 
494 // up/down of beam pipe
495                 halfx=xsize/2.;
496                 halfy=(ysize/2.-holeScaled/2.)/2.;          
497                 new TBRIK(nameSense,"Sens. region","void",halfx,halfy,0.25);
498                 
499                 node->cd();
500                 ypos=holeScaled/2.+((ysize/2.-holeScaled/2.)/2.);
501                 sprintf(nameNode,"MUON%d",200+id);
502                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense,0,ypos,0,"");
503                 nodeS->SetLineColor(kColorMUON2);
504                 
505                 node->cd();
506                 ypos=-1.*ypos;
507                 sprintf(nameNode,"MUON%d",300+id);
508                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense,0,ypos,0,"");
509                 nodeS->SetLineColor(kColorMUON2);
510                 
511 // left/right of beam pipe
512                 halfx=(xsize/2.-holeScaled/2.)/2.;
513                 halfy=holeScaled/2.;    
514                 new TBRIK(nameSense1,"Sens. region","void",halfx,halfy,0.25);
515                 
516                 node->cd();
517                 xpos=holeScaled/2.+((xsize/2.-holeScaled/2.)/2.);           
518                 sprintf(nameNode,"MUON%d",400+id);
519                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense1,xpos,0,0,"");
520                 nodeS->SetLineColor(kColorMUON2);
521                 
522                 node->cd();
523                 xpos=-1.*xpos;
524                 sprintf(nameNode,"MUON%d",500+id);
525                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense1,xpos,0,0,"");
526                 nodeS->SetLineColor(kColorMUON2);
527                 
528 // missing corners
529                 halfx=17.*zscale/2.;
530                 halfy=halfx;
531                 new TBRIK(nameSense2,"Sens. region","void",halfx,halfy,0.25);
532                 
533                 node->cd();
534                 xpos=holeScaled/2.-halfx;
535                 ypos=xpos;
536                 sprintf(nameNode,"MUON%d",600+id);
537                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
538                 nodeS->SetLineColor(kColorMUON2);
539                 
540                 node->cd();
541                 ypos=-1.*xpos;
542                 sprintf(nameNode,"MUON%d",700+id);
543                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
544                 nodeS->SetLineColor(kColorMUON2);
545                 
546                 node->cd();
547                 xpos=-1.*xpos;
548                 sprintf(nameNode,"MUON%d",800+id);
549                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
550                 nodeS->SetLineColor(kColorMUON2);
551                 
552                 node->cd();
553                 ypos=-1.*xpos;
554                 sprintf(nameNode,"MUON%d",900+id);
555                 nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
556                 nodeS->SetLineColor(kColorMUON2);
557             } 
558         }
559     }
560 }
561
562
563 //___________________________________________
564 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
565 {
566    return 9999;
567 }
568
569 //___________________________________________
570 void AliMUON::MakeBranch(Option_t* option)
571 {
572     // Create Tree branches for the MUON.
573     const Int_t kBufferSize = 4000;
574     char branchname[30];
575     sprintf(branchname,"%sCluster",GetName());
576     
577     AliDetector::MakeBranch(option);
578     
579     if (fPadHits   && gAlice->TreeH()) {
580         gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
581         printf("Making Branch %s for clusters\n",branchname);
582     }
583     
584 // one branch for digits per chamber
585     Int_t i;
586     
587     for (i=0; i<AliMUONConstants::NCh() ;i++) {
588         sprintf(branchname,"%sDigits%d",GetName(),i+1);
589         
590         if (fDchambers   && gAlice->TreeD()) {
591             gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
592             printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
593         }       
594     }
595     
596     printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
597     
598 // one branch for raw clusters per chamber
599     for (i=0; i<AliMUONConstants::NTrackingCh() ;i++) {
600         sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
601         
602         if (fRawClusters   && gAlice->TreeR()) {
603             gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
604             printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
605         }       
606     }
607     
608 // one branch for global trigger
609     sprintf(branchname,"%sGlobalTrigger",GetName());
610     if (fGlobalTrigger && gAlice->TreeR()) {  
611         gAlice->TreeR()->Branch(branchname,&fGlobalTrigger,kBufferSize);
612         printf("Making Branch %s for Global Trigger\n",branchname);
613     }
614 // one branch for local trigger
615     sprintf(branchname,"%sLocalTrigger",GetName());
616     if (fLocalTrigger && gAlice->TreeR()) {  
617         gAlice->TreeR()->Branch(branchname,&fLocalTrigger,kBufferSize);
618         printf("Making Branch %s for Local Trigger\n",branchname);
619     }
620     
621 }
622
623 //___________________________________________
624 void AliMUON::SetTreeAddress()
625 {
626   // Set branch address for the Hits and Digits Tree.
627   char branchname[30];
628   AliDetector::SetTreeAddress();
629
630   TBranch *branch;
631   TTree *treeH = gAlice->TreeH();
632   TTree *treeD = gAlice->TreeD();
633   TTree *treeR = gAlice->TreeR();
634
635   if (treeH) {
636     if (fPadHits) {
637       branch = treeH->GetBranch("MUONCluster");
638       if (branch) branch->SetAddress(&fPadHits);
639     }
640   }
641
642   if (treeD) {
643       for (int i=0; i<AliMUONConstants::NCh(); i++) {
644           sprintf(branchname,"%sDigits%d",GetName(),i+1);
645           if (fDchambers) {
646               branch = treeD->GetBranch(branchname);
647               if (branch) branch->SetAddress(&((*fDchambers)[i]));
648           }
649       }
650   }
651
652   // printf("SetTreeAddress --- treeR address  %p \n",treeR);
653
654   if (treeR) {
655       for (int i=0; i<AliMUONConstants::NTrackingCh(); i++) {
656           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
657           if (fRawClusters) {
658               branch = treeR->GetBranch(branchname);
659               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
660           }
661       }
662
663       if (fLocalTrigger) {
664         branch = treeR->GetBranch("MUONLocalTrigger");
665         if (branch) branch->SetAddress(&fLocalTrigger);
666       }
667       if (fGlobalTrigger) {
668         branch = treeR->GetBranch("MUONGlobalTrigger");
669         if (branch) branch->SetAddress(&fGlobalTrigger);
670       }
671   }
672 }
673 //___________________________________________
674 void AliMUON::ResetHits()
675 {
676   // Reset number of clusters and the cluster array for this detector
677   AliDetector::ResetHits();
678   fNPadHits = 0;
679   if (fPadHits) fPadHits->Clear();
680 }
681
682 //____________________________________________
683 void AliMUON::ResetDigits()
684 {
685     //
686     // Reset number of digits and the digits array for this detector
687     //
688     for ( int i=0;i<AliMUONConstants::NCh();i++ ) {
689         if ((*fDchambers)[i])    ((TClonesArray*)(*fDchambers)[i])->Clear();
690         if (fNdch)  fNdch[i]=0;
691     }
692 }
693 //____________________________________________
694 void AliMUON::ResetRawClusters()
695 {
696     //
697     // Reset number of raw clusters and the raw clust array for this detector
698     //
699     for ( int i=0;i<AliMUONConstants::NTrackingCh();i++ ) {
700         if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
701         if (fNrawch)  fNrawch[i]=0;
702     }
703 }
704
705 //____________________________________________
706 void AliMUON::ResetTrigger()
707 {
708   //  Reset Local and Global Trigger 
709   fNGlobalTrigger = 0;
710   if (fGlobalTrigger) fGlobalTrigger->Clear();
711   fNLocalTrigger = 0;
712   if (fLocalTrigger) fLocalTrigger->Clear();
713 }
714
715 //____________________________________________
716 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
717 {
718     Int_t i=2*(id-1);
719     ((AliMUONChamber*) (*fChambers)[i])  ->SetPadSize(isec,p1,p2);
720     ((AliMUONChamber*) (*fChambers)[i+1])->SetPadSize(isec,p1,p2);
721 }
722
723 //___________________________________________
724 void AliMUON::SetChambersZ(const Float_t *Z)
725 {
726   // Set Z values for all chambers (tracking and trigger)
727   // from the array pointed to by "Z"
728     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
729         ((AliMUONChamber*) ((*fChambers)[ch]))->SetZ(Z[ch]);
730     return;
731 }
732
733 //___________________________________________
734 void AliMUON::SetChambersZToDefault()
735 {
736   // Set Z values for all chambers (tracking and trigger)
737   // to default values
738   SetChambersZ(AliMUONConstants::DefaultChamberZ());
739   return;
740 }
741
742 //___________________________________________
743 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
744 {
745     Int_t i=2*(id-1);
746     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
747     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
748 }
749
750 //___________________________________________
751 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
752 {
753     Int_t i=2*(id-1);
754     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
755     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
756 }
757
758 //___________________________________________
759 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
760 {
761     Int_t i=2*(id-1);
762     ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
763     ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
764 }
765
766 //___________________________________________
767 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
768 {
769     Int_t i=2*(id-1);
770     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
771     ((AliMUONChamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
772 }
773
774 //___________________________________________
775 void AliMUON::SetMaxStepGas(Float_t p1)
776 {
777      fMaxStepGas=p1;
778 }
779
780 //___________________________________________
781 void AliMUON::SetMaxStepAlu(Float_t p1)
782 {
783     fMaxStepAlu=p1;
784 }
785
786 //___________________________________________
787 void AliMUON::SetMaxDestepGas(Float_t p1)
788 {
789     fMaxDestepGas=p1;
790 }
791
792 //___________________________________________
793 void AliMUON::SetMaxDestepAlu(Float_t p1)
794 {
795     fMaxDestepAlu=p1;
796 }
797 //___________________________________________
798 void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax)
799 {
800    fAccCut=acc;
801    fAccMin=angmin;
802    fAccMax=angmax;
803 }
804 //___________________________________________
805 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
806 {
807     ((AliMUONChamber*) (*fChambers)[id])->SetSegmentationModel(isec, segmentation);
808
809 }
810 //___________________________________________
811 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
812 {
813     ((AliMUONChamber*) (*fChambers)[id])->SetResponseModel(response);
814 }
815
816 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
817 {
818     ((AliMUONChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
819 }
820
821 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
822 {
823     ((AliMUONChamber*) (*fChambers)[id])->SetNsec(nsec);
824 }
825
826
827 //___________________________________________
828
829
830
831 void AliMUON::MakePadHits(Float_t xhit,Float_t yhit, Float_t zhit,
832                           Float_t eloss, Float_t tof,  Int_t idvol)
833 {
834 //
835 //  Calls the charge disintegration method of the current chamber and adds
836 //  the simulated cluster to the root treee 
837 //
838     Int_t clhits[7];
839     Float_t newclust[6][500];
840     Int_t nnew;
841     
842     
843 //
844 //  Integrated pulse height on chamber
845
846     
847     clhits[0]=fNhits+1;
848 //
849 //
850     ((AliMUONChamber*) (*fChambers)[idvol])
851         ->DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newclust);
852     Int_t ic=0;
853     
854 //
855 //  Add new clusters
856     for (Int_t i=0; i<nnew; i++) {
857         if (Int_t(newclust[3][i]) > 0) {
858             ic++;
859 // Cathode plane
860             clhits[1] = Int_t(newclust[5][i]);
861 //  Cluster Charge
862             clhits[2] = Int_t(newclust[0][i]);
863 //  Pad: ix
864             clhits[3] = Int_t(newclust[1][i]);
865 //  Pad: iy 
866             clhits[4] = Int_t(newclust[2][i]);
867 //  Pad: charge
868             clhits[5] = Int_t(newclust[3][i]);
869 //  Pad: chamber sector
870             clhits[6] = Int_t(newclust[4][i]);
871             
872             AddPadHit(clhits);
873         }
874     }
875 }
876
877 //----------------------------------------------------------------------
878
879 void AliMUON::Digitise(Int_t nev,Int_t bgrEvent,Option_t *option,Option_t *opt,Text_t *filename)
880 {
881     // keep galice.root for signal and name differently the file for 
882     // background when add! otherwise the track info for signal will be lost !
883   
884     static Bool_t first=kTRUE;
885     static TFile *file;
886     char *addBackground = strstr(option,"Add");
887
888
889     AliMUONChamber*   iChamber;
890     AliSegmentation*  segmentation;
891
892     
893     Int_t trk[50];
894     Int_t chtrk[50];  
895     TObjArray *list=new TObjArray;
896     static TClonesArray *pAddress=0;
897     if(!pAddress) pAddress=new TClonesArray("TVector",1000);
898     Int_t digits[5]; 
899
900     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
901     AliHitMap** hitMap= new AliHitMap* [AliMUONConstants::NCh()];
902     for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {hitMap[i]=0;}
903     if (addBackground ) {
904         if(first) {
905             fFileName=filename;
906             cout<<"filename"<<fFileName<<endl;
907             file=new TFile(fFileName);
908             cout<<"I have opened "<<fFileName<<" file "<<endl;
909             fHits2     = new TClonesArray("AliMUONHit",1000  );
910             fPadHits2 = new TClonesArray("AliMUONPadHit",10000);
911         }           
912         first=kFALSE;
913         file->cd();
914         //file->ls();
915         // Get Hits Tree header from file
916         if(fHits2) fHits2->Clear();
917         if(fPadHits2) fPadHits2->Clear();
918         if(fTrH1) delete fTrH1;
919         fTrH1=0;
920         
921         char treeName[20];
922         sprintf(treeName,"TreeH%d",bgrEvent);
923         fTrH1 = (TTree*)gDirectory->Get(treeName);
924         //printf("fTrH1 %p of treename %s for event %d \n",fTrH1,treeName,bgrEvent);
925         
926         if (!fTrH1) {
927             printf("ERROR: cannot find Hits Tree for event:%d\n",bgrEvent);
928         }
929         // Set branch addresses
930         TBranch *branch;
931         char branchname[20];
932         sprintf(branchname,"%s",GetName());
933         if (fTrH1 && fHits2) {
934             branch = fTrH1->GetBranch(branchname);
935             if (branch) branch->SetAddress(&fHits2);
936         }
937         if (fTrH1 && fPadHits2) {
938             branch = fTrH1->GetBranch("MUONCluster");
939             if (branch) branch->SetAddress(&fPadHits2);
940         }
941 // test
942         //Int_t ntracks1 =(Int_t)fTrH1->GetEntries();
943         //printf("background - ntracks1 - %d\n",ntracks1);
944     }
945     //
946     // loop over cathodes
947     //
948     AliHitMap* hm;
949     Int_t countadr=0;
950     for (int icat=0; icat<2; icat++) { 
951         Int_t counter=0;
952         Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
953         for (Int_t i =0; i<AliMUONConstants::NCh(); i++) {
954             iChamber=(AliMUONChamber*) (*fChambers)[i];
955             if (iChamber->Nsec()==1 && icat==1) {
956                 continue;
957             } else {
958                 segmentation=iChamber->SegmentationModel(icat+1);
959             }
960             hitMap[i] = new AliMUONHitMapA1(segmentation, list);
961             nmuon[i]=0;
962         }
963         //printf("Start loop over tracks \n");     
964 //
965 //   Loop over tracks
966 //
967
968         TTree *treeH = gAlice->TreeH();
969         Int_t ntracks =(Int_t) treeH->GetEntries();
970         Int_t jj;
971         Float_t ** xhit = new Float_t * [2];
972         for (jj=0; jj<2; jj++) xhit[jj] = new Float_t  [AliMUONConstants::NCh()];
973         Float_t ** yhit = new Float_t * [2];
974         for (jj=0; jj<2; jj++) yhit[jj] = new Float_t  [AliMUONConstants::NCh()];
975         
976         for (Int_t track=0; track<ntracks; track++) {
977             gAlice->ResetHits();
978             treeH->GetEvent(track);
979             
980 //
981 //   Loop over hits
982             for(AliMUONHit* mHit=(AliMUONHit*)pMUON->FirstHit(-1); 
983                 mHit;
984                 mHit=(AliMUONHit*)pMUON->NextHit()) 
985             {
986                 Int_t   nch   = mHit->fChamber-1;  // chamber number
987                 if (nch > AliMUONConstants::NCh()-1) continue;
988                 iChamber = &(pMUON->Chamber(nch));
989                 // new 17.07.99
990                 if (addBackground) {
991
992                     if (mHit->fParticle == kMuonPlus 
993                         || mHit->fParticle == kMuonMinus) {
994                         xhit[nch][nmuon[nch]]=mHit->fX;
995                         yhit[nch][nmuon[nch]]=mHit->fY;
996                         nmuon[nch]++;
997                         if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]);
998                     }
999                 }
1000                 
1001
1002
1003                 
1004 //
1005 // Loop over pad hits
1006                 for (AliMUONPadHit* mPad=
1007                          (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
1008                      mPad;
1009                      mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits))
1010                 {
1011                     Int_t cathode  = mPad->fCathode;    // cathode number
1012                     Int_t ipx      = mPad->fPadX;       // pad number on X
1013                     Int_t ipy      = mPad->fPadY;       // pad number on Y
1014                     Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
1015 //
1016 //
1017                     
1018                     if (cathode != (icat+1)) continue;
1019                     // fill the info array
1020                     Float_t thex, they, thez;
1021                     segmentation=iChamber->SegmentationModel(cathode);
1022                     segmentation->GetPadC(ipx,ipy,thex,they,thez);
1023 //                  Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1024 //                  if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1025
1026                     new((*pAddress)[countadr++]) TVector(2);
1027                     TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
1028                     trinfo(0)=(Float_t)track;
1029                     trinfo(1)=(Float_t)iqpad;
1030
1031                     digits[0]=ipx;
1032                     digits[1]=ipy;
1033                     digits[2]=iqpad;
1034                     digits[3]=iqpad;
1035                     if (mHit->fParticle == kMuonPlus ||
1036                         mHit->fParticle == kMuonMinus) {
1037                         digits[4]=mPad->fHitNumber;
1038                     } else digits[4]=-1;
1039
1040                     AliMUONTransientDigit* pdigit;
1041                     // build the list of fired pads and update the info
1042                     if (!hitMap[nch]->TestHit(ipx, ipy)) {
1043
1044                         list->AddAtAndExpand(
1045                             new AliMUONTransientDigit(nch,digits),counter);
1046                         
1047                         hitMap[nch]->SetHit(ipx, ipy, counter);
1048                         counter++;
1049                         pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
1050                         // list of tracks
1051                         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
1052                         trlist->Add(&trinfo);
1053                     } else {
1054                         pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
1055                         // update charge
1056                         (*pdigit).fSignal+=iqpad;
1057                         (*pdigit).fPhysics+=iqpad;                      
1058                         // update list of tracks
1059                         TObjArray* trlist=(TObjArray*)pdigit->TrackList();
1060                         Int_t lastEntry=trlist->GetLast();
1061                         TVector *pTrack=(TVector*)trlist->At(lastEntry);
1062                         TVector &ptrk=*pTrack;
1063                         Int_t lastTrack=Int_t(ptrk(0));
1064                         Int_t lastCharge=Int_t(ptrk(1));
1065                         if (lastTrack==track) {
1066                             lastCharge+=iqpad;
1067                             trlist->RemoveAt(lastEntry);
1068                             trinfo(0)=lastTrack;
1069                             trinfo(1)=lastCharge;
1070                             trlist->AddAt(&trinfo,lastEntry);
1071                         } else {
1072                             trlist->Add(&trinfo);
1073                         }
1074                         // check the track list
1075                         Int_t nptracks=trlist->GetEntriesFast();
1076                         if (nptracks > 2) {
1077                             for (Int_t tr=0;tr<nptracks;tr++) {
1078                                 TVector *ppTrack=(TVector*)trlist->At(tr);
1079                                 TVector &pptrk=*ppTrack;
1080                                 trk[tr]=Int_t(pptrk(0));
1081                                 chtrk[tr]=Int_t(pptrk(1));
1082                             }
1083                         } // end if nptracks
1084                     } //  end if pdigit
1085                 } //end loop over clusters
1086             } // hit loop
1087         } // track loop
1088
1089         // open the file with background
1090        
1091         if (addBackground) {
1092             ntracks =(Int_t)fTrH1->GetEntries();
1093 //
1094 //   Loop over tracks
1095 //
1096             for (Int_t track=0; track<ntracks; track++) {
1097
1098                 if (fHits2)       fHits2->Clear();
1099                 if (fPadHits2)   fPadHits2->Clear();
1100
1101                 fTrH1->GetEvent(track);
1102 //
1103 //   Loop over hits
1104                 AliMUONHit* mHit;
1105                 for(int i=0;i<fHits2->GetEntriesFast();++i) 
1106                 {       
1107                     mHit=(AliMUONHit*) (*fHits2)[i];
1108                     Int_t   nch   = mHit->fChamber-1;  // chamber number
1109                     if (nch >9) continue;
1110                     iChamber = &(pMUON->Chamber(nch));
1111                     Int_t rmin = (Int_t)iChamber->RInner();
1112                     Int_t rmax = (Int_t)iChamber->ROuter();
1113                     Float_t xbgr=mHit->fX;
1114                     Float_t ybgr=mHit->fY;
1115                     Bool_t cond=kFALSE;
1116                     
1117                     for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
1118                         Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
1119                             +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
1120                         if (dist<100) cond=kTRUE;
1121                     }
1122                     if (!cond) continue;
1123                     
1124 //
1125 // Loop over pad hits
1126                     for (AliMUONPadHit* mPad=
1127                              (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits2);
1128                          mPad;
1129                          mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits2))
1130                     {
1131                         //                  mPad = (AliMUONPadHit*) (*fPadHits2)[j];
1132                         Int_t cathode  = mPad->fCathode;    // cathode number
1133                         Int_t ipx      = mPad->fPadX;       // pad number on X
1134                         Int_t ipy      = mPad->fPadY;       // pad number on Y
1135                         Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
1136
1137                         if (cathode != (icat+1)) continue;
1138                         Float_t thex, they, thez;
1139                         segmentation=iChamber->SegmentationModel(cathode);
1140                         segmentation->GetPadC(ipx,ipy,thex,they,thez);
1141                         Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1142                         if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1143                         new((*pAddress)[countadr++]) TVector(2);
1144                         TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
1145                         trinfo(0)=-1;  // tag background
1146                         trinfo(1)=-1;
1147                         
1148                         digits[0]=ipx;
1149                         digits[1]=ipy;
1150                         digits[2]=iqpad;
1151                         digits[3]=0;
1152                         digits[4]=-1;
1153                         
1154                         AliMUONTransientDigit* pdigit;
1155                         // build the list of fired pads and update the info
1156                         if (!hitMap[nch]->TestHit(ipx, ipy)) {
1157                             list->AddAtAndExpand(new AliMUONTransientDigit(nch,digits),counter);
1158                             
1159                             hitMap[nch]->SetHit(ipx, ipy, counter);
1160                             counter++;
1161                             
1162                             pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
1163                             // list of tracks
1164                             TObjArray *trlist=(TObjArray*)pdigit->
1165                                 TrackList();
1166                             trlist->Add(&trinfo);
1167                         } else {
1168                             pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
1169                             // update charge
1170                             (*pdigit).fSignal+=iqpad;
1171                             
1172                             // update list of tracks
1173                             TObjArray* trlist=(TObjArray*)pdigit->
1174                                 TrackList();
1175                             Int_t lastEntry=trlist->GetLast();
1176                             TVector *pTrack=(TVector*)trlist->
1177                                 At(lastEntry);
1178                             TVector &ptrk=*pTrack;
1179                             Int_t lastTrack=Int_t(ptrk(0));
1180                             if (lastTrack==-1) {
1181                                 continue;
1182                             } else {
1183                                 trlist->Add(&trinfo);
1184                             }
1185                                 // check the track list
1186                             Int_t nptracks=trlist->GetEntriesFast();
1187                             if (nptracks > 0) {
1188                                 for (Int_t tr=0;tr<nptracks;tr++) {
1189                                     TVector *ppTrack=(TVector*)trlist->At(tr);
1190                                     TVector &pptrk=*ppTrack;
1191                                     trk[tr]=Int_t(pptrk(0));
1192                                     chtrk[tr]=Int_t(pptrk(1));
1193                                 }
1194                             } // end if nptracks
1195                         } //  end if pdigit
1196                     } //end loop over clusters
1197                 } // hit loop
1198             } // track loop
1199             //Int_t nentr2=list->GetEntriesFast();
1200             //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1201             TTree *fAli=gAlice->TreeK();
1202             TFile *file=NULL;
1203             
1204             if (fAli) file =fAli->GetCurrentFile();
1205             file->cd();
1206         } // if addBackground
1207         delete [] xhit;
1208         delete [] yhit;
1209
1210         Int_t tracks[10];
1211         Int_t charges[10];
1212         Int_t nentries=list->GetEntriesFast();
1213         
1214         for (Int_t nent=0;nent<nentries;nent++) {
1215             AliMUONTransientDigit *address=(AliMUONTransientDigit*)list->At(nent);
1216             if (address==0) continue; 
1217             Int_t ich=address->fChamber;
1218             Int_t q=address->fSignal; 
1219             iChamber=(AliMUONChamber*) (*fChambers)[ich];
1220 //
1221 //  Digit Response (noise, threshold, saturation, ...)
1222 //              if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise; 
1223             AliMUONResponse * response=iChamber->ResponseModel();
1224             q=response->DigitResponse(q);
1225             
1226             if (!q) continue;
1227             
1228             digits[0]=address->fPadX;
1229             digits[1]=address->fPadY;
1230             digits[2]=q;
1231             digits[3]=address->fPhysics;
1232             digits[4]=address->fHit;
1233             
1234             TObjArray* trlist=(TObjArray*)address->TrackList();
1235             Int_t nptracks=trlist->GetEntriesFast();
1236             //printf("nptracks, trlist   %d  %p\n",nptracks,trlist);
1237
1238             // this was changed to accomodate the real number of tracks
1239             if (nptracks > 10) {
1240                 cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1241                 nptracks=10;
1242             }
1243             if (nptracks > 2) {
1244                 printf("Attention - nptracks > 2  %d \n",nptracks);
1245                 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1246             }
1247             for (Int_t tr=0;tr<nptracks;tr++) {
1248                 TVector *ppP=(TVector*)trlist->At(tr);
1249                 if(!ppP ) printf("ppP - %p\n",ppP);
1250                 TVector &pp  =*ppP;
1251                 tracks[tr]=Int_t(pp(0));
1252                 charges[tr]=Int_t(pp(1));
1253                 //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
1254             }      //end loop over list of tracks for one pad
1255             // Sort list of tracks according to charge
1256             if (nptracks > 1) {
1257                 SortTracks(tracks,charges,nptracks);
1258             }
1259             if (nptracks < 10 ) {
1260                 for (Int_t i=nptracks; i<10; i++) {
1261                     tracks[i]=0;
1262                     charges[i]=0;
1263                 }
1264             }
1265             
1266             // fill digits
1267             pMUON->AddDigits(ich,tracks,charges,digits);
1268             // delete trlist;
1269         }
1270         //cout<<"I'm out of the loops for digitisation"<<endl;
1271         //      gAlice->GetEvent(nev);
1272         gAlice->TreeD()->Fill();
1273         pMUON->ResetDigits();
1274         list->Delete();
1275
1276         
1277         for(Int_t ii=0;ii<AliMUONConstants::NCh();++ii) {
1278             if (hitMap[ii]) {
1279                 hm=hitMap[ii];
1280                 delete hm;
1281                 hitMap[ii]=0;
1282             }
1283         }
1284         delete [] nmuon;    
1285     } //end loop over cathodes
1286     delete [] hitMap;
1287     char hname[30];
1288     sprintf(hname,"TreeD%d",nev);
1289     gAlice->TreeD()->Write(hname);
1290     // reset tree
1291     gAlice->TreeD()->Reset();
1292     delete list;
1293     
1294     pAddress->Delete();
1295     // gObjectTable->Print();
1296 }
1297
1298 void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
1299 {
1300   //
1301   // Sort the list of tracks contributing to a given digit
1302   // Only the 3 most significant tracks are acctually sorted
1303   //
1304   
1305   //
1306   //  Loop over signals, only 3 times
1307   //
1308   
1309   Int_t qmax;
1310   Int_t jmax;
1311   Int_t idx[3] = {-2,-2,-2};
1312   Int_t jch[3] = {-2,-2,-2};
1313   Int_t jtr[3] = {-2,-2,-2};
1314   Int_t i,j,imax;
1315   
1316   if (ntr<3) imax=ntr;
1317   else imax=3;
1318   for(i=0;i<imax;i++){
1319     qmax=0;
1320     jmax=0;
1321     
1322     for(j=0;j<ntr;j++){
1323       
1324       if((i == 1 && j == idx[i-1]) 
1325          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
1326       
1327       if(charges[j] > qmax) {
1328         qmax = charges[j];
1329         jmax=j;
1330       }       
1331     } 
1332     
1333     if(qmax > 0) {
1334       idx[i]=jmax;
1335       jch[i]=charges[jmax]; 
1336       jtr[i]=tracks[jmax]; 
1337     }
1338     
1339   } 
1340   
1341   for(i=0;i<3;i++){
1342     if (jtr[i] == -2) {
1343          charges[i]=0;
1344          tracks[i]=0;
1345     } else {
1346          charges[i]=jch[i];
1347          tracks[i]=jtr[i];
1348     }
1349   }
1350
1351 }
1352
1353 //___________________________________________
1354 void AliMUON::Trigger(Int_t nev){
1355 // call the Trigger Algorithm and fill TreeR
1356
1357   Int_t singlePlus[3]  = {0,0,0}; 
1358   Int_t singleMinus[3] = {0,0,0}; 
1359   Int_t singleUndef[3] = {0,0,0};
1360   Int_t pairUnlike[3]  = {0,0,0}; 
1361   Int_t pairLike[3]    = {0,0,0};
1362
1363   ResetTrigger();
1364   AliMUONTriggerDecision* decision= new AliMUONTriggerDecision(1);
1365   decision->Trigger();   
1366   decision->GetGlobalTrigger(singlePlus, singleMinus, singleUndef,
1367                              pairUnlike, pairLike);
1368 // add a local trigger in the list 
1369   AddGlobalTrigger(singlePlus, singleMinus, singleUndef, pairUnlike, pairLike);
1370   Int_t i;
1371   
1372   for (Int_t icirc=0; icirc<AliMUONConstants::NTriggerCircuit(); icirc++) { 
1373       if(decision->GetITrigger(icirc)==1) {
1374           Int_t localtr[7]={0,0,0,0,0,0,0};      
1375           Int_t loLpt[2]={0,0}; Int_t loHpt[2]={0,0}; Int_t loApt[2]={0,0};
1376           decision->GetLutOutput(icirc, loLpt, loHpt, loApt);
1377           localtr[0] = icirc;
1378           localtr[1] = decision->GetStripX11(icirc);
1379           localtr[2] = decision->GetDev(icirc);
1380           localtr[3] = decision->GetStripY11(icirc);
1381           for (i=0; i<2; i++) {    // convert the Lut output in 1 digit 
1382               localtr[4] = localtr[4]+Int_t(loLpt[i]*TMath::Power(2,i));
1383               localtr[5] = localtr[5]+Int_t(loHpt[i]*TMath::Power(2,i));
1384               localtr[6] = localtr[6]+Int_t(loApt[i]*TMath::Power(2,i));
1385           }
1386           AddLocalTrigger(localtr);  // add a local trigger in the list
1387       }
1388   }
1389   delete decision;
1390
1391   gAlice->TreeR()->Fill();
1392   ResetTrigger();
1393   char hname[30];
1394   sprintf(hname,"TreeR%d",nev);
1395   gAlice->TreeR()->Write(hname);
1396   gAlice->TreeR()->Reset();
1397   printf("\n End of trigger for event %d", nev);
1398 }
1399
1400
1401 //____________________________________________
1402 void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
1403 {
1404     TClonesArray *dig1, *dig2;
1405     Int_t ndig, k;
1406     dig1 = new TClonesArray("AliMUONDigit",1000);
1407     dig2 = new TClonesArray("AliMUONDigit",1000);
1408     AliMUONDigit *digit;
1409 //
1410 // Loop on chambers and on cathode planes
1411 //
1412     
1413     for (Int_t ich=0;ich<10;ich++) {
1414         AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
1415         AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();    
1416         gAlice->ResetDigits();
1417         gAlice->TreeD()->GetEvent(lastEntry);
1418         TClonesArray *muonDigits  = this->DigitsAddress(ich);
1419         ndig=muonDigits->GetEntriesFast();
1420         printf("\n 1 Found %d digits in %p %d", ndig, muonDigits,ich);
1421         TClonesArray &lhits1 = *dig1;
1422         Int_t n=0;
1423         for (k=0; k<ndig; k++) {
1424             digit=(AliMUONDigit*) muonDigits->UncheckedAt(k);
1425             if (rec->TestTrack(digit->fTracks[0]))
1426                 new(lhits1[n++]) AliMUONDigit(*digit);
1427         }
1428         gAlice->ResetDigits();
1429         gAlice->TreeD()->GetEvent(lastEntry+1);
1430         muonDigits  = this->DigitsAddress(ich);
1431         ndig=muonDigits->GetEntriesFast();
1432         printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
1433         TClonesArray &lhits2 = *dig2;
1434         n=0;
1435         
1436         for (k=0; k<ndig; k++) {
1437             digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
1438             if (rec->TestTrack(digit->fTracks[0]))
1439             new(lhits2[n++]) AliMUONDigit(*digit);
1440         }
1441
1442         if (rec) {       
1443             AliMUONClusterInput::Instance()->SetDigits(ich, dig1, dig2);
1444             rec->FindRawClusters();
1445         }
1446         dig1->Delete();
1447         dig2->Delete();
1448     } // for ich
1449     gAlice->TreeR()->Fill();
1450     ResetRawClusters();
1451     char hname[30];
1452     sprintf(hname,"TreeR%d",nev);
1453     gAlice->TreeR()->Write(hname);
1454     gAlice->TreeR()->Reset();
1455     printf("\n End of cluster finding for event %d", nev);
1456     
1457     delete dig1;
1458     delete dig2;
1459     //gObjectTable->Print();
1460 }
1461  
1462
1463 void AliMUON::Streamer(TBuffer &R__b)
1464 {
1465    // Stream an object of class AliMUON.
1466       AliMUONChamber        *iChamber;
1467       AliMUONTriggerCircuit *iTriggerCircuit;
1468       AliSegmentation       *segmentation;
1469       AliMUONResponse       *response;
1470       TClonesArray          *digitsaddress;
1471       TClonesArray          *rawcladdress;
1472       Int_t i;
1473       if (R__b.IsReading()) {
1474           Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1475           AliDetector::Streamer(R__b);
1476           R__b >> fNPadHits;
1477           R__b >> fPadHits; // diff
1478           R__b >> fNLocalTrigger;       
1479           R__b >> fLocalTrigger;       
1480           R__b >> fNGlobalTrigger;       
1481           R__b >> fGlobalTrigger;   
1482           R__b >> fDchambers;
1483           R__b >> fRawClusters;
1484           R__b.ReadArray(fNdch);
1485           R__b.ReadArray(fNrawch);
1486           R__b >> fAccCut;
1487           R__b >> fAccMin;
1488           R__b >> fAccMax; 
1489           R__b >> fChambers;
1490           R__b >> fTriggerCircuits;
1491           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
1492               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
1493               iTriggerCircuit->Streamer(R__b);
1494           }
1495 // Stream chamber related information
1496           for (i =0; i<AliMUONConstants::NCh(); i++) {
1497               iChamber=(AliMUONChamber*) (*fChambers)[i];
1498               iChamber->Streamer(R__b);
1499               if (iChamber->Nsec()==1) {
1500                   segmentation=iChamber->SegmentationModel(1);
1501                   if (segmentation)
1502                       segmentation->Streamer(R__b);
1503               } else {
1504                   segmentation=iChamber->SegmentationModel(1);
1505                   if (segmentation)
1506                       segmentation->Streamer(R__b);
1507                   if (segmentation)
1508                       segmentation=iChamber->SegmentationModel(2);
1509                   segmentation->Streamer(R__b);
1510               }
1511               response=iChamber->ResponseModel();
1512               if (response)
1513                   response->Streamer(R__b);       
1514               digitsaddress=(TClonesArray*) (*fDchambers)[i];
1515               digitsaddress->Streamer(R__b);
1516               if (i < AliMUONConstants::NTrackingCh()) {
1517                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1518                   rawcladdress->Streamer(R__b);
1519               }
1520           }
1521           
1522       } else {
1523           R__b.WriteVersion(AliMUON::IsA());
1524           AliDetector::Streamer(R__b);
1525           R__b << fNPadHits;
1526           R__b << fPadHits; // diff
1527           R__b << fNLocalTrigger;       
1528           R__b << fLocalTrigger;       
1529           R__b << fNGlobalTrigger;       
1530           R__b << fGlobalTrigger; 
1531           R__b << fDchambers;
1532           R__b << fRawClusters;
1533           R__b.WriteArray(fNdch, AliMUONConstants::NCh());
1534           R__b.WriteArray(fNrawch, AliMUONConstants::NTrackingCh());
1535           
1536           R__b << fAccCut;
1537           R__b << fAccMin;
1538           R__b << fAccMax; 
1539           
1540           R__b << fChambers;
1541           R__b << fTriggerCircuits;
1542           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
1543               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
1544               iTriggerCircuit->Streamer(R__b);
1545           }
1546           for (i =0; i<AliMUONConstants::NCh(); i++) {
1547               iChamber=(AliMUONChamber*) (*fChambers)[i];
1548               iChamber->Streamer(R__b);
1549               if (iChamber->Nsec()==1) {
1550                   segmentation=iChamber->SegmentationModel(1);
1551                   if (segmentation)
1552                       segmentation->Streamer(R__b);
1553               } else {
1554                   segmentation=iChamber->SegmentationModel(1);
1555                   if (segmentation)
1556                       segmentation->Streamer(R__b);
1557                   segmentation=iChamber->SegmentationModel(2);
1558                   if (segmentation)
1559                       segmentation->Streamer(R__b);
1560               }
1561               response=iChamber->ResponseModel();
1562               if (response)
1563                   response->Streamer(R__b);
1564               digitsaddress=(TClonesArray*) (*fDchambers)[i];
1565               digitsaddress->Streamer(R__b);
1566               if (i < AliMUONConstants::NTrackingCh()) {
1567                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1568                   rawcladdress->Streamer(R__b);
1569               }
1570           }
1571       }
1572 }
1573 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters) 
1574 {
1575 //
1576     // Initialise the pad iterator
1577     // Return the address of the first padhit for hit
1578     TClonesArray *theClusters = clusters;
1579     Int_t nclust = theClusters->GetEntriesFast();
1580     if (nclust && hit->fPHlast > 0) {
1581         AliMUON::fMaxIterPad=hit->fPHlast;
1582         AliMUON::fCurIterPad=hit->fPHfirst;
1583         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
1584     } else {
1585         return 0;
1586     }
1587 }
1588
1589 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters) 
1590 {
1591     AliMUON::fCurIterPad++;
1592     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
1593         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
1594     } else {
1595         return 0;
1596     }
1597 }
1598
1599
1600 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
1601 {
1602     TClonesArray *muonRawCluster  = RawClustAddress(ichamber);
1603     ResetRawClusters();
1604     TTree *treeR = gAlice->TreeR();
1605     Int_t nent=(Int_t)treeR->GetEntries();
1606     treeR->GetEvent(nent-2+icathod-1);
1607     //treeR->GetEvent(icathod);
1608     //Int_t nrawcl = (Int_t)muonRawCluster->GetEntriesFast();
1609
1610     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)muonRawCluster->UncheckedAt(icluster);
1611     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
1612     
1613     return  mRaw;
1614 }
1615
1616 AliMUON& AliMUON::operator = (const AliMUON& rhs)
1617 {
1618 // copy operator
1619 // dummy version
1620     return *this;
1621 }
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640