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