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