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