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