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