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