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