MUON/AliMUON.cxx

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