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