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 /* $Id$ */
  17
  18
  19 ///////////////////////////////////////////////
  20 //  Manager and hits classes for set:MUON     //
  21 ////////////////////////////////////////////////
  22
  23 #include "Riostream.h"
  24
  25 #include <AliPDG.h>
  26 #include <TBRIK.h>
  27 #include <TCanvas.h>
  28 #include <TDirectory.h>
  29 #include <TFile.h>
  30 #include <TGeometry.h>
  31 #include <TMinuit.h>
  32 #include <TNode.h>
  33 #include <TNtuple.h>
  34 #include <TObjArray.h>
  35 #include <TObject.h>
  36 #include <TObjectTable.h>
  37 #include <TPad.h>
  38 #include <TParticle.h>
  39 #include <TROOT.h>
  40 #include <TRandom.h>
  41 #include <TRotMatrix.h>
  42 #include <TTUBE.h>
  43 #include <TTUBE.h>
  44 #include <TTree.h>
  45 #include <TVector.h>
  46 #include <TVirtualMC.h>
  47
  48 #include "AliConst.h"
  49 #include "AliHeader.h"
  50 #include "AliHitMap.h"
  51 #include "AliLoader.h"
  52 #include "AliRunDigitizer.h"
  53 #include "AliESD.h"
  54 #include "AliESDMuonTrack.h"
  55 #include "AliMC.h"
  56 #include "AliMUONLoader.h"
  57 #include "AliMUON.h"
  58 #include "AliMUONTriggerTrack.h"
  59 #include "AliMUONEventReconstructor.h"
  60 #include "AliMUONClusterReconstructor.h"
  61 #include "AliMUONTrack.h"
  62 #include "AliMUONTrackParam.h"
  63 #include "AliMUONChamberTrigger.h"
  64 #include "AliMUONClusterFinderAZ.h"
  65 #include "AliMUONClusterInput.h"
  66 #include "AliMUONConstants.h"
  67 #include "AliMUONDigit.h"
  68 #include "AliMUONGlobalTrigger.h"
  69 #include "AliMUONHit.h"
  70 #include "AliMUONHitMapA1.h"
  71 #include "AliMUONLocalTrigger.h"
  72 #include "AliMUONMerger.h"
  73 #include "AliMUONPadHit.h"
  74 #include "AliMUONRawCluster.h"
  75 #include "AliMUONTransientDigit.h"
  76 #include "AliMUONTriggerCircuit.h"
  77 #include "AliMUONTriggerDecision.h"
  78 #include "AliMUONVGeometryBuilder.h"
  79 #include "AliRun.h"
  80 #include "AliMUONDigitizerv2.h"
  81 #include "AliMUONSDigitizerv1.h"
  82
  83
  84 // Defaults parameters for Z positions of chambers
  85 // taken from values for "stations" in AliMUON::AliMUON
  86 //     const Float_t zch[7]={528, 690., 975., 1249., 1449., 1610, 1710.};
  87 // and from array "dstation" in AliMUONv1::CreateGeometry
  88 //          Float_t dstation[5]={20., 20., 20, 20., 20.};
  89 //     for tracking chambers,
  90 //          according to (Z1 = zch - dstation) and  (Z2 = zch + dstation)
  91 //          for the first and second chambers in the station, respectively,
  92 // and from "DTPLANES" in AliMUONv1::CreateGeometry
  93 //           const Float_t DTPLANES = 15.;
  94 //     for trigger chambers,
  95 //          according to (Z1 = zch) and  (Z2 = zch + DTPLANES)
  96 //          for the first and second chambers in the station, respectively
  97
  98 ClassImp(AliMUON)
  99 //__________________________________________________________________
 100 AliMUON::AliMUON()
 101 {
 102 // Default Constructor
 103 //
 104     fNCh             = 0;
 105     fNTrackingCh     = 0;
 106     fIshunt          = 0;
 107     fChambers        = 0;
 108     fGeometryBuilders = 0;
 109     fTriggerCircuits = 0;
 110     fAccMin          = 0.;
 111     fAccMax          = 0.;
 112     fAccCut          = kFALSE;
 113     fMerger          = 0;
 114     fMUONData        = 0;
 115     fSplitLevel      = 0;
 116 }
 117 //__________________________________________________________________
 118 AliMUON::AliMUON(const char *name, const char *title)
 119   : AliDetector(name,title)
 120 {
 121 //Begin_Html
 122 /*
 123 <img src="gif/alimuon.gif">
 124 */
 125 //End_Html
 126   fMUONData  = 0x0;
 127   fSplitLevel= 0;
 128   fIshunt     =  0;
 129
 130   fNCh             = AliMUONConstants::NCh();
 131   fNTrackingCh     = AliMUONConstants::NTrackingCh();
 132
 133   SetMarkerColor(kRed);//
 134 //
 135 // Creating List of Chambers
 136     Int_t ch;
 137     fChambers = new TObjArray(AliMUONConstants::NCh());
 138     fGeometryBuilders = new TObjArray(AliMUONConstants::NCh());
 139
 140     // Loop over stations
 141     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
 142       // Loop over 2 chambers in the station
 143       for (Int_t stCH = 0; stCH < 2; stCH++) {
 144         //
 145         //
 146         //    Default Parameters for Muon Tracking Stations
 147         ch = 2 * st + stCH;
 148         if (ch < AliMUONConstants::NTrackingCh()) {
 149           fChambers->AddAt(new AliMUONChamber(ch),ch);
 150         } else {
 151           fChambers->AddAt(new AliMUONChamberTrigger(ch),ch);
 152         }
 153         AliMUONChamber* chamber = (AliMUONChamber*) fChambers->At(ch);
 154         //chamber->SetGid(0);
 155         // Default values for Z of chambers
 156         chamber->SetZ(AliMUONConstants::DefaultChamberZ(ch));
 157         //
 158         chamber->InitGeo(AliMUONConstants::DefaultChamberZ(ch));
 159         //          Set chamber inner and outer radius to default
 160         chamber->SetRInner(AliMUONConstants::Dmin(st)/2);
 161         chamber->SetROuter(AliMUONConstants::Dmax(st)/2);
 162         //
 163       } // Chamber stCH (0, 1) in
 164     }     // Station st (0...)
 165
 166     // Negatives values are ignored by geant3 CONS200 in the calculation of the tracking parameters
 167     fMaxStepGas=0.1;
 168     fMaxStepAlu=0.1;
 169     fMaxDestepGas=-1;
 170     fMaxDestepAlu=-1;
 171
 172     fMaxIterPad   = 0;
 173     fCurIterPad   = 0;
 174
 175     fAccMin          = 0.;
 176     fAccMax          = 0.;
 177     fAccCut          = kFALSE;
 178
 179     // cp new design of AliMUONTriggerDecision
 180     fTriggerCircuits = new TObjArray(AliMUONConstants::NTriggerCircuit());
 181     for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
 182       fTriggerCircuits->AddAt(new AliMUONTriggerCircuit(),circ);
 183     }
 184     fMerger = 0;
 185 }
 186 //____________________________________________________________________
 187 AliMUON::AliMUON(const AliMUON& rMUON):AliDetector(rMUON)
 188 {
 189 // Dummy copy constructor
 190     ;
 191
 192 }
 193 //____________________________________________________________________
 194 AliMUON::~AliMUON()
 195 {
 196 // Destructor
 197   if(fDebug) printf("%s: Calling AliMUON destructor !!!\n",ClassName());
 198   fIshunt  = 0;
 199   if (fMerger) delete fMerger;
 200
 201   if (fGeometryBuilders){
 202     fGeometryBuilders->Delete();
 203     delete fGeometryBuilders;
 204   }
 205   if (fChambers){
 206     fChambers->Delete();
 207     delete fChambers;
 208   }
 209   if (fTriggerCircuits){
 210     fTriggerCircuits->Delete();
 211     delete fTriggerCircuits;
 212   }
 213   delete fMUONData;
 214 }
 215 //_____________________________________________________________________________
 216 void AliMUON::AddGeometryBuilder(AliMUONVGeometryBuilder* geomBuilder)
 217 {
 218 // Adds the geometry builder to the list
 219 // ---
 220
 221   fGeometryBuilders->Add(geomBuilder);
 222 }
 223 //____________________________________________________________________
 224 void AliMUON::BuildGeometry()
 225 {
 226 // Geometry for event display
 227   for (Int_t i=0; i<7; i++) {
 228     for (Int_t j=0; j<2; j++) {
 229       Int_t id=2*i+j+1;
 230       this->Chamber(id-1).SegmentationModel(1)->Draw("eventdisplay");
 231     }
 232   }
 233 }
 234 //___________________________________________________________________
 235 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
 236 {
 237   return 9999;
 238 }
 239 //__________________________________________________________________
 240 void  AliMUON::SetTreeAddress()
 241 {
 242   GetMUONData()->SetLoader(fLoader);
 243   //  GetMUONData()->MakeBranch("D,S,RC");
 244   //  GetMUONData()->SetTreeAddress("H,D,S,RC");
 245   GetMUONData()->SetTreeAddress("H");
 246   if (fHits !=  GetMUONData()->Hits())  {
 247     if ( gAlice->GetMCApp() )
 248       if ( gAlice->GetMCApp()->GetHitLists() ) {
 249         fHits = GetMUONData()->Hits();
 250         gAlice->GetMCApp()->AddHitList(fHits); // For purifyKine, only necessary when Hit list is created in AliMUONData
 251       }
 252   }
 253   fHits = GetMUONData()->Hits(); // Added by Ivana to use the methods FisrtHit, NextHit of AliDetector
 254 }
 255
 256 //____________________________________________________________________
 257 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 258 {
 259 // Set the pad size for chamber id and cathode isec
 260     Int_t i=2*(id-1);
 261     ((AliMUONChamber*) fChambers->At(i))  ->SetPadSize(isec,p1,p2);
 262     ((AliMUONChamber*) fChambers->At(i+1))->SetPadSize(isec,p1,p2);
 263 }
 264
 265 //___________________________________________
 266 void AliMUON::SetChambersZ(const Float_t *Z)
 267 {
 268   // Set Z values for all chambers (tracking and trigger)
 269   // from the array pointed to by "Z"
 270     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
 271         ((AliMUONChamber*) fChambers->At(ch))->SetZ(Z[ch]);
 272     return;
 273 }
 274 //_________________________________________________________________
 275 void AliMUON::SetChambersZToDefault()
 276 {
 277   // Set Z values for all chambers (tracking and trigger)
 278   // to default values
 279   SetChambersZ(AliMUONConstants::DefaultChamberZ());
 280   return;
 281 }
 282 //_________________________________________________________________
 283 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
 284 {
 285 // Set the inverse charge slope for chamber id
 286     Int_t i=2*(id-1);    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
 287     //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
 288     ((AliMUONChamber*) fChambers->At(i))->SetChargeSlope(p1);
 289     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSlope(p1);
 290 }
 291 //__________________________________________________________________
 292 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 293 {
 294 // Set sigma of charge spread for chamber id
 295     Int_t i=2*(id-1);
 296     ((AliMUONChamber*) fChambers->At(i))->SetChargeSpread(p1,p2);
 297     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSpread(p1,p2);
 298 }
 299 //___________________________________________________________________
 300 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
 301 {
 302 // Set integration limits for charge spread
 303     Int_t i=2*(id-1);
 304     ((AliMUONChamber*) fChambers->At(i))->SetSigmaIntegration(p1);
 305     ((AliMUONChamber*) fChambers->At(i+1))->SetSigmaIntegration(p1);
 306 }
 307
 308 //__________________________________________________________________
 309 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
 310 {
 311 // Set maximum number for ADCcounts (saturation)
 312     Int_t i=2*(id-1);
 313     ((AliMUONChamber*) fChambers->At(i))->SetMaxAdc(p1);
 314     ((AliMUONChamber*) fChambers->At(i+1))->SetMaxAdc(p1);
 315 }
 316 //__________________________________________________________________
 317 void AliMUON::SetMaxStepGas(Float_t p1)
 318 {
 319 // Set stepsize in gas
 320   fMaxStepGas=p1;
 321 }
 322 //__________________________________________________________________
 323 void AliMUON::SetMaxStepAlu(Float_t p1)
 324 {
 325 // Set step size in Alu
 326     fMaxStepAlu=p1;
 327 }
 328 //__________________________________________________________________
 329 void AliMUON::SetMaxDestepGas(Float_t p1)
 330 {
 331 // Set maximum step size in Gas
 332     fMaxDestepGas=p1;
 333 }
 334 //__________________________________________________________________
 335 void AliMUON::SetMaxDestepAlu(Float_t p1)
 336 {
 337 // Set maximum step size in Alu
 338   fMaxDestepAlu=p1;
 339 }
 340 //___________________________________________________________________
 341 void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
 342 {
 343 // Set acceptance cuts
 344   fAccCut=acc;
 345   fAccMin=angmin*TMath::Pi()/180;
 346   fAccMax=angmax*TMath::Pi()/180;
 347   Int_t ch;
 348   if (acc) {
 349     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
 350       // Loop over 2 chambers in the station
 351       for (Int_t stCH = 0; stCH < 2; stCH++) {
 352         ch = 2 * st + stCH;
 353         //         Set chamber inner and outer radius according to acceptance cuts
 354         Chamber(ch).SetRInner(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMin)));
 355         Chamber(ch).SetROuter(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMax)));
 356       } // chamber loop
 357     } // station loop
 358   }
 359 }
 360
 361 //____________________________________________________________________
 362 Float_t  AliMUON::GetMaxStepGas() const
 363 {
 364 // Return stepsize in gas
 365
 366   return fMaxStepGas;
 367 }
 368
 369 //____________________________________________________________________
 370 Float_t  AliMUON::GetMaxStepAlu() const
 371 {
 372 // Return step size in Alu
 373
 374   return fMaxStepAlu;
 375 }
 376
 377 //____________________________________________________________________
 378 Float_t  AliMUON::GetMaxDestepGas() const
 379 {
 380 // Return maximum step size in Gas
 381
 382   return fMaxDestepGas;
 383 }
 384
 385 //____________________________________________________________________
 386 Float_t  AliMUON::GetMaxDestepAlu() const
 387 {
 388 // Return maximum step size in Gas
 389
 390   return fMaxDestepAlu;
 391 }
 392
 393 //____________________________________________________________________
 394 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
 395 {
 396 // Set the segmentation for chamber id cathode isec
 397     ((AliMUONChamber*) fChambers->At(id))->SetSegmentationModel(isec, segmentation);
 398
 399 }
 400 //____________________________________________________________________
 401 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
 402 {
 403 // Set the response for chamber id
 404     ((AliMUONChamber*) fChambers->At(id))->SetResponseModel(response);
 405 }
 406 //____________________________________________________________________
 407 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
 408 {
 409 // Set ClusterFinder for chamber id
 410     ((AliMUONChamber*) fChambers->At(id))->SetReconstructionModel(reconst);
 411 }
 412 //____________________________________________________________________
 413 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
 414 {
 415 // Set number of segmented cathods for chamber id
 416     ((AliMUONChamber*) fChambers->At(id))->SetNsec(nsec);
 417 }
 418 //____________________________________________________________________
 419 AliDigitizer* AliMUON::CreateDigitizer(AliRunDigitizer* manager) const
 420 {
 421   return new AliMUONDigitizerv2(manager);
 422 }
 423 //_____________________________________________________________________
 424 void AliMUON::SDigits2Digits()
 425 {
 426
 427 // write TreeD here
 428
 429     if (!fMerger) {
 430       if (gAlice->GetDebug()>0) {
 431         cerr<<"AliMUON::SDigits2Digits: create default AliMUONMerger "<<endl;
 432         cerr<<" no merging, just digitization of 1 event will be done"<<endl;
 433       }
 434       fMerger = new AliMUONMerger();
 435     }
 436     fMerger->Init();
 437     fMerger->Digitise();
 438     char hname[30];
 439     //    sprintf(hname,"TreeD%d",fLoader->GetHeader()->GetEvent());
 440     fLoader->TreeD()->Write(hname,TObject::kOverwrite);
 441     fLoader->TreeD()->Reset();
 442 }
 443
 444 //_____________________________________________________________________
 445 void AliMUON::Hits2SDigits()
 446 {
 447   // Adaption of AliMUONSDigitizerv1 to be excuted by the AliSimulation framework
 448   AliRunLoader* runLoader = fLoader->GetRunLoader();
 449   AliRunDigitizer   * manager = new AliRunDigitizer(1,1);
 450   manager->SetInputStream(0,runLoader->GetFileName(),AliConfig::GetDefaultEventFolderName());
 451   AliMUONDigitizer * dMUON   = new AliMUONSDigitizerv1(manager);
 452   fLoader->LoadHits("READ");
 453   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
 454     runLoader->GetEvent(iEvent);
 455     dMUON->Exec("");
 456   }
 457   fLoader->UnloadHits();
 458 }
 459 //_______________________________________________________________________
 460 AliLoader* AliMUON::MakeLoader(const char* topfoldername)
 461 {
 462 //builds standard getter (AliLoader type)
 463 //if detector wants to use castomized getter, it must overload this method
 464
 465  if (GetDebug())
 466    Info("MakeLoader",
 467         "Creating standard getter for detector %s. Top folder is %s.",
 468          GetName(),topfoldername);
 469  fLoader   = new AliLoader(GetName(),topfoldername);
 470  fMUONData = new AliMUONData(fLoader,GetName(),GetName());
 471  fMUONData->SetSplitLevel(fSplitLevel);
 472  return fLoader;
 473 }
 474
 475 //_______________________________________________________________________
 476 void AliMUON::Trigger(Int_t /*nev*/){
 477 // call the Trigger Algorithm and fill TreeR
 478
 479   Int_t singlePlus[3]  = {0,0,0};
 480   Int_t singleMinus[3] = {0,0,0};
 481   Int_t singleUndef[3] = {0,0,0};
 482   Int_t pairUnlike[3]  = {0,0,0};
 483   Int_t pairLike[3]    = {0,0,0};
 484
 485   ResetTrigger();
 486   AliMUONTriggerDecision* decision= new AliMUONTriggerDecision(fLoader,1);
 487   decision->Trigger();
 488   decision->GetGlobalTrigger(singlePlus, singleMinus, singleUndef,
 489                              pairUnlike, pairLike);
 490
 491   // add a local trigger in the list
 492   GetMUONData()->AddGlobalTrigger(singlePlus, singleMinus, singleUndef, pairUnlike, pairLike);
 493   Int_t i;
 494
 495   for (Int_t icirc=0; icirc<AliMUONConstants::NTriggerCircuit(); icirc++) {
 496     if(decision->GetITrigger(icirc)==1) {
 497       Int_t localtr[7]={0,0,0,0,0,0,0};
 498       Int_t loLpt[2]={0,0}; Int_t loHpt[2]={0,0}; Int_t loApt[2]={0,0};
 499       decision->GetLutOutput(icirc, loLpt, loHpt, loApt);
 500       localtr[0] = icirc;
 501       localtr[1] = decision->GetStripX11(icirc);
 502       localtr[2] = decision->GetDev(icirc);
 503       localtr[3] = decision->GetStripY11(icirc);
 504       for (i=0; i<2; i++) {    // convert the Lut output in 1 digit
 505         localtr[4] = localtr[4]+Int_t(loLpt[i]*TMath::Power(2,i));
 506         localtr[5] = localtr[5]+Int_t(loHpt[i]*TMath::Power(2,i));
 507         localtr[6] = localtr[6]+Int_t(loApt[i]*TMath::Power(2,i));
 508       }
 509       GetMUONData()->AddLocalTrigger(localtr);  // add a local trigger in the list
 510     }
 511   }
 512
 513   delete decision;
 514
 515   //  fLoader->TreeR()->Fill();
 516   GetMUONData()->Fill("GLT"); //Filling Global and Local Trigger GLT
 517   //  char hname[30];
 518   //  sprintf(hname,"TreeR%d",nev);
 519   //  fLoader->TreeR()->Write(hname,TObject::kOverwrite);
 520     //  fLoader->TreeR()->Reset();
 521   fLoader->WriteRecPoints("OVERWRITE");
 522
 523   //  printf("\n End of trigger for event %d\n", nev);
 524 }
 525
 526 //____________________________________________________________________
 527 void AliMUON::Digits2Reco()
 528 {
 529   FindClusters();
 530   Int_t nev = gAlice->GetHeader()->GetEvent();
 531   GetMUONData()->Fill("RC"); //Filling Reconstructed Cluster
 532   fLoader->WriteRecPoints("OVERWRITE");
 533   GetMUONData()->ResetRawClusters();
 534   Info("Digits2Reco","End of cluster finding for event %d", nev);
 535 }
 536 //____________________________________________________________________
 537 void AliMUON::FindClusters()
 538 {
 539 //
 540 //  Perform cluster finding
 541 //
 542     TClonesArray *dig1, *dig2;
 543     Int_t ndig, k;
 544     dig1 = new TClonesArray("AliMUONDigit",1000);
 545     dig2 = new TClonesArray("AliMUONDigit",1000);
 546     AliMUONDigit *digit;
 547 // Loop on chambers and on cathode planes
 548 //
 549     ResetRawClusters();
 550     TClonesArray * muonDigits;
 551
 552     for (Int_t ich = 0; ich < 10; ich++) {
 553       //PH      AliMUONChamber* iChamber = (AliMUONChamber*) (*fChambers)[ich];
 554         AliMUONChamber* iChamber = (AliMUONChamber*) fChambers->At(ich);
 555         AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();
 556
 557         ResetDigits();
 558         GetMUONData()->GetCathode(0);
 559         //TClonesArray *
 560         muonDigits = GetMUONData()->Digits(ich);
 561         ndig=muonDigits->GetEntriesFast();
 562         if(fDebug)
 563         printf("\n 1 Found %d digits in %p chamber %d", ndig, muonDigits,ich);
 564         TClonesArray &lhits1 = *dig1;
 565         Int_t n = 0;
 566         for (k = 0; k < ndig; k++) {
 567             digit = (AliMUONDigit*) muonDigits->UncheckedAt(k);
 568             if (rec->TestTrack(digit->Track(0)))
 569                 new(lhits1[n++]) AliMUONDigit(*digit);
 570         }
 571         GetMUONData()->ResetDigits();
 572         GetMUONData()->GetCathode(1);
 573         muonDigits =  GetMUONData()->Digits(ich);
 574         ndig=muonDigits->GetEntriesFast();
 575         if(fDebug)
 576         printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
 577         TClonesArray &lhits2 = *dig2;
 578         n=0;
 579
 580         for (k=0; k<ndig; k++) {
 581             digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
 582             if (rec->TestTrack(digit->Track(0)))
 583             new(lhits2[n++]) AliMUONDigit(*digit);
 584         }
 585
 586         if (rec) {
 587             AliMUONClusterInput::Instance()->SetDigits(ich, dig1, dig2);
 588             rec->FindRawClusters();
 589         }
 590         dig1->Delete();
 591         dig2->Delete();
 592     } // for ich
 593     delete dig1;
 594     delete dig2;
 595 }
 596 //______________________________________________________________________
 597 #ifdef never
 598 void AliMUON::Streamer(TBuffer &R__b)_
 599 {
 600    // Stream an object of class AliMUON.
 601       AliMUONChamber        *iChamber;
 602       AliMUONTriggerCircuit *iTriggerCircuit;
 603       AliSegmentation       *segmentation;
 604       AliMUONResponse       *response;
 605       TClonesArray          *digitsaddress;
 606       TClonesArray          *rawcladdress;
 607       Int_t i;
 608       if (R__b.IsReading()) {
 609           Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 610           AliDetector::Streamer(R__b);
 611           R__b >> fNPadHits;
 612           R__b >> fPadHits; // diff
 613           R__b >> fNLocalTrigger;
 614           R__b >> fLocalTrigger;
 615           R__b >> fNGlobalTrigger;
 616           R__b >> fGlobalTrigger;
 617           R__b >> fDchambers;
 618           R__b >> fRawClusters;
 619           R__b.ReadArray(fNdch);
 620           R__b.ReadArray(fNrawch);
 621           R__b >> fAccCut;
 622           R__b >> fAccMin;
 623           R__b >> fAccMax;
 624           R__b >> fChambers;
 625           R__b >> fTriggerCircuits;
 626           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
 627               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
 628               iTriggerCircuit->Streamer(R__b);
 629           }
 630 // Stream chamber related information
 631           for (i =0; i<AliMUONConstants::NCh(); i++) {
 632               iChamber=(AliMUONChamber*) (*fChambers)[i];
 633               iChamber->Streamer(R__b);
 634               if (iChamber->Nsec()==1) {
 635                   segmentation=iChamber->SegmentationModel(1);
 636                   if (segmentation)
 637                       segmentation->Streamer(R__b);
 638               } else {
 639                   segmentation=iChamber->SegmentationModel(1);
 640                   if (segmentation)
 641                       segmentation->Streamer(R__b);
 642                   if (segmentation)
 643                       segmentation=iChamber->SegmentationModel(2);
 644                   segmentation->Streamer(R__b);
 645               }
 646               response=iChamber->ResponseModel();
 647               if (response)
 648                   response->Streamer(R__b);
 649               digitsaddress=(TClonesArray*) (*fDchambers)[i];
 650               digitsaddress->Streamer(R__b);
 651               if (i < AliMUONConstants::NTrackingCh()) {
 652                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
 653                   rawcladdress->Streamer(R__b);
 654               }
 655           }
 656
 657       } else {
 658           R__b.WriteVersion(AliMUON::IsA());
 659           AliDetector::Streamer(R__b);
 660           R__b << fNPadHits;
 661           R__b << fPadHits; // diff
 662           R__b << fNLocalTrigger;
 663           R__b << fLocalTrigger;
 664           R__b << fNGlobalTrigger;
 665           R__b << fGlobalTrigger;
 666           R__b << fDchambers;
 667           R__b << fRawClusters;
 668           R__b.WriteArray(fNdch, AliMUONConstants::NCh());
 669           R__b.WriteArray(fNrawch, AliMUONConstants::NTrackingCh());
 670
 671           R__b << fAccCut;
 672           R__b << fAccMin;
 673           R__b << fAccMax;
 674
 675           R__b << fChambers;
 676           R__b << fTriggerCircuits;
 677           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
 678               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
 679               iTriggerCircuit->Streamer(R__b);
 680           }
 681           for (i =0; i<AliMUONConstants::NCh(); i++) {
 682               iChamber=(AliMUONChamber*) (*fChambers)[i];
 683               iChamber->Streamer(R__b);
 684               if (iChamber->Nsec()==1) {
 685                   segmentation=iChamber->SegmentationModel(1);
 686                   if (segmentation)
 687                       segmentation->Streamer(R__b);
 688               } else {
 689                   segmentation=iChamber->SegmentationModel(1);
 690                   if (segmentation)
 691                       segmentation->Streamer(R__b);
 692                   segmentation=iChamber->SegmentationModel(2);
 693                   if (segmentation)
 694                       segmentation->Streamer(R__b);
 695               }
 696               response=iChamber->ResponseModel();
 697               if (response)
 698                   response->Streamer(R__b);
 699               digitsaddress=(TClonesArray*) (*fDchambers)[i];
 700               digitsaddress->Streamer(R__b);
 701               if (i < AliMUONConstants::NTrackingCh()) {
 702                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
 703                   rawcladdress->Streamer(R__b);
 704               }
 705           }
 706       }
 707 }
 708 #endif
 709 //_______________________________________________________________________
 710 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters)
 711 {
 712 // to be removed
 713     // Initialise the pad iterator
 714     // Return the address of the first padhit for hit
 715     TClonesArray *theClusters = clusters;
 716     Int_t nclust = theClusters->GetEntriesFast();
 717     if (nclust && hit->PHlast() > 0) {
 718         AliMUON::fMaxIterPad=hit->PHlast();
 719         AliMUON::fCurIterPad=hit->PHfirst();
 720         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
 721     } else {
 722         return 0;
 723     }
 724 }
 725 //_______________________________________________________________________
 726 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters)
 727 {
 728   // To be removed
 729 // Get next pad (in iterator)
 730 //
 731     AliMUON::fCurIterPad++;
 732     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
 733         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
 734     } else {
 735         return 0;
 736     }
 737 }
 738 //_______________________________________________________________________
 739
 740 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
 741 {
 742 //
 743 //  Return rawcluster (icluster) for chamber ichamber and cathode icathod
 744 //  Obsolete ??
 745     TClonesArray *muonRawCluster  = GetMUONData()->RawClusters(ichamber);
 746     ResetRawClusters();
 747     TTree *treeR = fLoader->TreeR();
 748     Int_t nent=(Int_t)treeR->GetEntries();
 749     treeR->GetEvent(nent-2+icathod-1);
 750     //treeR->GetEvent(icathod);
 751     //Int_t nrawcl = (Int_t)muonRawCluster->GetEntriesFast();
 752
 753     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)muonRawCluster->UncheckedAt(icluster);
 754     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
 755
 756     return  mRaw;
 757 }
 758 //________________________________________________________________________
 759 void   AliMUON::SetMerger(AliMUONMerger* merger)
 760 {
 761 // Set pointer to merger
 762     fMerger = merger;
 763 }
 764 //________________________________________________________________________
 765 AliMUONMerger*  AliMUON::Merger()
 766 {
 767 // Return pointer to merger
 768     return fMerger;
 769 }
 770 //________________________________________________________________________
 771 AliMUON& AliMUON::operator = (const AliMUON& /*rhs*/)
 772 {
 773 // copy operator
 774 // dummy version
 775     return *this;
 776 }
 777 //________________________________________________________________________
 778 void AliMUON::Reconstruct() const
 779 {
 780
 781 //  AliLoader* loader = GetLoader();
 782
 783   AliRunLoader* runLoader = fLoader->GetRunLoader();
 784   Int_t nEvents = runLoader->GetNumberOfEvents();
 785
 786 // used local container for each method
 787 // passing fLoader as argument, could be avoided ???
 788   AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(fLoader);
 789   AliMUONData* dataEvent = recoEvent->GetMUONData();
 790
 791   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(fLoader);
 792   AliMUONData* dataCluster = recoCluster->GetMUONData();
 793
 794   AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(fLoader);
 795   AliMUONData* dataTrig = trigDec->GetMUONData();
 796
 797
 798   for (Int_t i = 0; i < 10; i++) {
 799     AliMUONClusterFinderVS *RecModel = new AliMUONClusterFinderVS();
 800     RecModel->SetGhostChi2Cut(10);
 801     recoCluster->SetReconstructionModel(i,RecModel);
 802   }
 803
 804   fLoader->LoadDigits("READ");
 805   fLoader->LoadRecPoints("RECREATE");
 806   fLoader->LoadTracks("RECREATE");
 807
 808   //   Loop over events
 809   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
 810     printf("Event %d\n",ievent);
 811     runLoader->GetEvent(ievent);
 812
 813     //----------------------- digit2cluster & Digits2Trigger -------------------
 814     if (!fLoader->TreeR()) fLoader->MakeRecPointsContainer();
 815
 816     // tracking branch
 817     dataCluster->MakeBranch("RC");
 818     dataCluster->SetTreeAddress("D,RC");
 819     recoCluster->Digits2Clusters();
 820     dataCluster->Fill("RC");
 821
 822     // trigger branch
 823     dataTrig->MakeBranch("GLT");
 824     dataTrig->SetTreeAddress("D,GLT");
 825     trigDec->Digits2Trigger();
 826     dataTrig->Fill("GLT");
 827
 828     fLoader->WriteRecPoints("OVERWRITE");
 829
 830     //---------------------------- Track & TriggerTrack ---------------------
 831     if (!fLoader->TreeT()) fLoader->MakeTracksContainer();
 832
 833     // trigger branch
 834     dataEvent->MakeBranch("RL"); //trigger track
 835     dataEvent->SetTreeAddress("RL");
 836     recoEvent->EventReconstructTrigger();
 837     dataEvent->Fill("RL");
 838
 839     // tracking branch
 840     dataEvent->MakeBranch("RT"); //track
 841     dataEvent->SetTreeAddress("RT");
 842     recoEvent->EventReconstruct();
 843     dataEvent->Fill("RT");
 844
 845     fLoader->WriteTracks("OVERWRITE");
 846
 847     //--------------------------- Resetting branches -----------------------
 848     dataCluster->ResetDigits();
 849     dataCluster->ResetRawClusters();
 850
 851     dataTrig->ResetDigits();
 852     dataTrig->ResetTrigger();
 853
 854     dataEvent->ResetRawClusters();
 855     dataEvent->ResetTrigger();
 856     dataEvent->ResetRecTracks();
 857     dataEvent->ResetRecTriggerTracks();
 858
 859   }
 860   fLoader->UnloadDigits();
 861   fLoader->UnloadRecPoints();
 862   fLoader->UnloadTracks();
 863
 864   delete recoCluster;
 865   delete recoEvent;
 866   delete trigDec;
 867 }
 868 //________________________________________________________________________
 869 void AliMUON::FillESD(AliESD* event) const
 870 {
 871
 872   TClonesArray* recTracksArray = 0;
 873   TClonesArray* recTrigTracksArray = 0;
 874
 875   //YS AliLoader* loader = GetLoader();
 876   AliRunLoader* runLoader = fLoader->GetRunLoader();
 877   fLoader->LoadTracks("READ"); //YS
 878
 879
 880   // declaration
 881   Int_t iEvent;
 882   Int_t nTrackHits;
 883   Double_t fitFmin;
 884
 885
 886   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
 887   Double_t xRec, yRec, zRec, chi2MatchTrigger;
 888   Bool_t matchTrigger;
 889
 890   //YS Int_t nEvents = runLoader->GetNumberOfEvents();
 891
 892   // setting pointer for tracks, triggertracks& trackparam at vertex
 893   AliMUONTrack* recTrack = 0;
 894   AliMUONTrackParam* trackParam = 0;
 895   AliMUONTriggerTrack* recTriggerTrack = 0;
 896
 897   iEvent = runLoader->GetEventNumber() ; //YS, seems not to be implemented yet (Ch. F)
 898   runLoader->GetEvent(iEvent);
 899
 900   // setting ESD MUON class
 901   AliESDMuonTrack* ESDTrack = new  AliESDMuonTrack() ;
 902
 903   //-------------------- trigger tracks-------------
 904   Long_t trigPat = 0;
 905   fMUONData->SetTreeAddress("RL");
 906   fMUONData->GetRecTriggerTracks();
 907   recTrigTracksArray = fMUONData->RecTriggerTracks();
 908
 909   // ready global trigger pattern from first track
 910   if (recTrigTracksArray)
 911     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
 912   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
 913
 914   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 915
 916   // -------------------- tracks-------------
 917   fMUONData->SetTreeAddress("RT");
 918   fMUONData->GetRecTracks();
 919   recTracksArray = fMUONData->RecTracks();
 920
 921   Int_t nRecTracks = 0;
 922   if (recTracksArray)
 923     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
 924
 925   // loop over tracks
 926   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 927
 928     // reading info from tracks
 929     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 930
 931     trackParam = recTrack->GetTrackParamAtVertex();
 932
 933     bendingSlope            = trackParam->GetBendingSlope();
 934     nonBendingSlope         = trackParam->GetNonBendingSlope();
 935     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
 936     xRec  = trackParam->GetNonBendingCoor();
 937     yRec  = trackParam->GetBendingCoor();
 938     zRec  = trackParam->GetZ();
 939
 940     nTrackHits       = recTrack->GetNTrackHits();
 941     fitFmin          = recTrack->GetFitFMin();
 942     matchTrigger     = recTrack->GetMatchTrigger();
 943     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 944
 945     // setting data member of ESD MUON
 946     ESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
 947     ESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
 948     ESDTrack->SetThetaY(TMath::ATan(bendingSlope));
 949     ESDTrack->SetZ(zRec);
 950     ESDTrack->SetBendingCoor(yRec);
 951     ESDTrack->SetNonBendingCoor(xRec);
 952     ESDTrack->SetChi2(fitFmin);
 953     ESDTrack->SetNHit(nTrackHits);
 954     ESDTrack->SetMatchTrigger(matchTrigger);
 955     ESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 956
 957     // storing ESD MUON Track into ESD Event
 958     if (nRecTracks != 0)
 959       event->AddMuonTrack(ESDTrack);
 960   } // end loop tracks
 961
 962   // add global trigger pattern
 963   if (nRecTracks != 0)
 964     event->SetTrigger(trigPat);
 965
 966   // reset muondata
 967   fMUONData->ResetRecTracks();
 968   fMUONData->ResetRecTriggerTracks();
 969
 970   //} // end loop on event
 971   fLoader->UnloadTracks();
 972 }
 973