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