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