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