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