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