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