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