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