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