MUON/AliMUONDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2000, 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 #include <Riostream.h>
  19 #include <TDirectory.h>
  20 #include <TPDGCode.h>
  21
  22 #include "AliRunLoader.h"
  23 #include "AliLoader.h"
  24 #include <Riostream.h>
  25
  26 #include "AliMUONDigitizer.h"
  27 #include "AliMUONConstants.h"
  28 #include "AliMUONChamber.h"
  29 #include "AliMUONHitMapA1.h"
  30 #include "AliMUON.h"
  31 #include "AliMUONChamber.h"
  32 #include "AliMUONConstants.h"
  33 #include "AliMUONDigit.h"
  34 #include "AliMUONDigitizer.h"
  35 #include "AliMUONHit.h"
  36 #include "AliMUONHitMapA1.h"
  37 #include "AliMUONPadHit.h"
  38 #include "AliMUONTransientDigit.h"
  39 #include "AliRun.h"
  40 #include "AliRunDigitizer.h"
  41
  42 ClassImp(AliMUONDigitizer)
  43
  44 //___________________________________________
  45 AliMUONDigitizer::AliMUONDigitizer() :AliDigitizer()
  46 {
  47 // Default ctor - don't use it
  48   fHits = 0;
  49   fPadHits = 0;
  50   fHitMap = 0;
  51   fTDList = 0;
  52 }
  53
  54 //___________________________________________
  55 AliMUONDigitizer::AliMUONDigitizer(AliRunDigitizer* manager)
  56     :AliDigitizer(manager)
  57 {
  58 // ctor which should be used
  59   fHitMap  = 0;
  60   fTDList  = 0;
  61   fHits    = 0;
  62   fPadHits = 0;
  63   fDebug   = 0;
  64   if (GetDebug()>2)
  65     cerr<<"AliMUONDigitizer::AliMUONDigitizer"
  66        <<"(AliRunDigitizer* manager) was processed"<<endl;
  67 }
  68
  69 //------------------------------------------------------------------------
  70 AliMUONDigitizer::~AliMUONDigitizer()
  71 {
  72 // Destructor
  73   if (fHits)       delete fHits;
  74   if (fPadHits)    delete fPadHits;
  75 }
  76
  77 //------------------------------------------------------------------------
  78 Bool_t AliMUONDigitizer::Exists(const AliMUONPadHit *padhit) const
  79 {
  80   return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
  81 }
  82
  83 //------------------------------------------------------------------------
  84 void AliMUONDigitizer::Update(AliMUONPadHit *padhit)
  85 {
  86     AliMUONTransientDigit *pdigit =
  87       static_cast<AliMUONTransientDigit*>(
  88       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
  89
  90     // update charge
  91     //
  92     Int_t iqpad    = padhit->QPad();        // charge per pad
  93     pdigit->AddSignal(iqpad);
  94     pdigit->AddPhysicsSignal(iqpad);
  95
  96     // update list of tracks
  97     //
  98     Int_t track, charge;
  99     track = fTrack+fMask;
 100     if (fSignal) {
 101       charge = iqpad;
 102     } else {
 103       charge = kBgTag;
 104     }
 105     pdigit->UpdateTrackList(track,charge);
 106 }
 107
 108 //------------------------------------------------------------------------
 109 void AliMUONDigitizer::CreateNew(AliMUONPadHit *padhit)
 110 {
 111 // Create new AliMUONTransientDigit and add it to the fTDList
 112
 113   fTDList->AddAtAndExpand(
 114     new AliMUONTransientDigit(fNch,fDigits),fCounter);
 115   fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
 116   AliMUONTransientDigit* pdigit =
 117     (AliMUONTransientDigit*)fTDList->Last();
 118   // list of tracks
 119   Int_t track, charge;
 120   track = fTrack+fMask;
 121   if (fSignal) {
 122     charge = padhit->QPad();
 123   } else {
 124     charge = kBgTag;
 125   }
 126   pdigit->AddToTrackList(track,charge);
 127   fCounter++;
 128 }
 129
 130
 131 //------------------------------------------------------------------------
 132 Bool_t AliMUONDigitizer::Init()
 133 {
 134 // Initialization
 135
 136   fHits     = new TClonesArray("AliMUONHit",1000);
 137   fPadHits  = new TClonesArray("AliMUONPadHit",1000);
 138   return kTRUE;
 139 }
 140
 141
 142 //------------------------------------------------------------------------
 143 //void AliMUONDigitizer::Digitize()
 144 void AliMUONDigitizer::Exec(Option_t* option)
 145 {
 146
 147   TString optionString = option;
 148   if (optionString.Data() == "deb") {
 149     cout<<"AliMUONDigitizer::Exec: called with option deb "<<endl;
 150     fDebug = 3;
 151   }
 152   AliMUONChamber*   iChamber;
 153   AliSegmentation*  segmentation;
 154
 155   if (GetDebug()>2) cerr<<"   AliMUONDigitizer::Digitize() starts"<<endl;
 156   fTDList = new TObjArray;
 157
 158   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 159   if (!pMUON) {
 160     cerr<<"AliMUONDigitizer::Digitize Error:"
 161        <<" module MUON not found in the input file"<<endl;
 162     return;
 163   }
 164
 165   AliRunLoader *rl, *orl;
 166   AliLoader *gime, *ogime;
 167   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 168   ogime = orl->GetLoader("MUONLoader");
 169
 170   pMUON->MakeBranchInTreeD(ogime->TreeD());
 171
 172   fHitMap= new AliMUONHitMapA1* [AliMUONConstants::NCh()];
 173
 174   //
 175   // loop over cathodes
 176   //
 177
 178   for (int icat = 0; icat < 2; icat++) {
 179     fCounter = 0;
 180     for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
 181       iChamber = &(pMUON->Chamber(i));
 182 //      if (!(iChamber->Nsec() == 1 && icat == 1)) {
 183        segmentation = iChamber->SegmentationModel(icat+1);
 184        fHitMap[i] = new AliMUONHitMapA1(segmentation, fTDList);
 185 //      }
 186     }
 187
 188
 189 // Loop over files to digitize
 190     fSignal = kTRUE;
 191     for (Int_t inputFile=0; inputFile<fManager->GetNinputs();
 192         inputFile++) {
 193 // Connect MUON branches
 194
 195       if (inputFile > 0 ) fSignal = kFALSE;
 196       TBranch *branchHits = 0;
 197       TBranch *branchPadHits = 0;
 198
 199       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
 200       gime = rl->GetLoader("MUONLoader");
 201
 202
 203       TTree *treeH = gime->TreeH();
 204       if (GetDebug()>2) {
 205        cerr<<"   inputFile , cathode = "<<inputFile<<" "
 206            <<icat<<endl;
 207        cerr<<"   treeH, fHits "<<treeH<<" "<<fHits<<endl;
 208       }
 209       if (treeH && fHits) {
 210        branchHits = treeH->GetBranch("MUON");
 211        if (branchHits) {
 212          fHits->Clear();
 213          branchHits->SetAddress(&fHits);
 214        }
 215        else
 216          Error("Exec","branch MUON was not found");
 217       }
 218       if (GetDebug()>2) cerr<<"   branchHits = "<<branchHits<<endl;
 219
 220       if (treeH && fPadHits) {
 221        branchPadHits = treeH->GetBranch("MUONCluster");
 222        if (branchPadHits)
 223          branchPadHits->SetAddress(&fPadHits);
 224        else
 225          Error("Exec","branch MUONCluster was not found");
 226       }
 227       if (GetDebug()>2) cerr<<"   branchPadHits = "<<branchPadHits<<endl;
 228
 229 //
 230 //   Loop over tracks
 231 //
 232
 233       Int_t ntracks = (Int_t) treeH->GetEntries();
 234
 235       for (fTrack = 0; fTrack < ntracks; fTrack++) {
 236        if (GetDebug()>2) cerr<<"   fTrack = "<<fTrack<<endl;
 237        fHits->Clear();
 238        fPadHits->Clear();
 239        branchHits->GetEntry(fTrack);
 240        branchPadHits->GetEntry(fTrack);
 241
 242 //
 243 //   Loop over hits
 244
 245        AliMUONHit* mHit;
 246        for(Int_t i = 0; i < fHits->GetEntriesFast(); ++i) {
 247          mHit = static_cast<AliMUONHit*>(fHits->At(i));
 248          fNch = mHit->Chamber()-1;  // chamber number
 249          if (fNch > AliMUONConstants::NCh()-1) {
 250            cerr<<"AliMUONDigitizer: ERROR: "
 251               <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
 252               <<fNch<<", "<< AliMUONConstants::NCh()<<endl;
 253            return;
 254          }
 255          iChamber = &(pMUON->Chamber(fNch));
 256 //
 257 // Loop over pad hits
 258          for (AliMUONPadHit* mPad =
 259                (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
 260               mPad;
 261               mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHits))
 262          {
 263            Int_t cathode  = mPad->Cathode();      // cathode number
 264            Int_t ipx      = mPad->PadX();         // pad number on X
 265            Int_t ipy      = mPad->PadY();         // pad number on Y
 266            Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
 267            if (cathode != (icat+1)) continue;
 268
 269            fMask = fManager->GetMask(inputFile);
 270            fDigits[0] = ipx;
 271            fDigits[1] = ipy;
 272            fDigits[2] = icat;
 273            fDigits[3] = iqpad;
 274            if (inputFile == 0) {
 275              fDigits[4] = iqpad;
 276            } else {
 277              fDigits[4] = 0;
 278            }
 279            if (mHit->Particle() == kMuonPlus ||
 280               mHit->Particle() == kMuonMinus) {
 281              fDigits[5] = (mPad->HitNumber()) + fMask;
 282            } else fDigits[5] = -1;
 283
 284            // build the list of fired pads and update the info,
 285            // fDigits is input for Update(mPad)
 286
 287            if (!Exists(mPad)) {
 288              CreateNew(mPad);
 289            } else {
 290              Update(mPad);
 291            } //  end if Exists(mPad)
 292          } //end loop over clusters
 293        } // hit loop
 294       } // track loop
 295     } // end file loop
 296     if (GetDebug()>2) cerr<<"END OF FILE LOOP"<<endl;
 297
 298     Int_t tracks[kMAXTRACKS];
 299     Int_t charges[kMAXTRACKS];
 300     Int_t nentries = fTDList->GetEntriesFast();
 301
 302     for (Int_t nent = 0; nent < nentries; nent++) {
 303       AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
 304       if (address == 0) continue;
 305       Int_t ich = address->Chamber();
 306       Int_t   q = address->Signal();
 307       iChamber = &(pMUON->Chamber(ich));
 308 //
 309 //  Digit Response (noise, threshold, saturation, ...)
 310       AliMUONResponse * response = iChamber->ResponseModel();
 311       q = response->DigitResponse(q,address);
 312
 313       if (!q) continue;
 314
 315       fDigits[0] = address->PadX();
 316       fDigits[1] = address->PadY();
 317       fDigits[2] = address->Cathode();
 318       fDigits[3] = q;
 319       fDigits[4] = address->Physics();
 320       fDigits[5] = address->Hit();
 321
 322       Int_t nptracks = address->GetNTracks();
 323
 324       if (nptracks > kMAXTRACKS) {
 325        if (GetDebug() >0) {
 326          cerr<<"AliMUONDigitizer: nptracks > 10 "<<nptracks;
 327          cerr<<"reset to max value "<<kMAXTRACKS<<endl;
 328        }
 329        nptracks = kMAXTRACKS;
 330       }
 331       if (nptracks > 2 && GetDebug() >2) {
 332        cerr<<"AliMUONDigitizer: nptracks > 2 "<<nptracks<<endl;
 333        printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
 334       }
 335       for (Int_t tr = 0; tr < nptracks; tr++) {
 336        tracks[tr]   = address->GetTrack(tr);
 337        charges[tr]  = address->GetCharge(tr);
 338       }      //end loop over list of tracks for one pad
 339       // Sort list of tracks according to charge
 340       if (nptracks > 1) {
 341        SortTracks(tracks,charges,nptracks);
 342       }
 343       if (nptracks < kMAXTRACKS ) {
 344        for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
 345          tracks[i]  = 0;
 346          charges[i] = 0;
 347        }
 348       }
 349
 350       // fill digits
 351       pMUON->AddDigits(ich,tracks,charges,fDigits);
 352     }
 353
 354 //    fManager->GetTreeD()->Fill();
 355     ogime->TreeD()->Fill();
 356
 357
 358     pMUON->ResetDigits();  //
 359     fTDList->Clear();
 360
 361
 362     for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
 363       if (fHitMap[ii]) {
 364        delete fHitMap[ii];
 365        fHitMap[ii] = 0;
 366       }
 367     }
 368   } //end loop over cathodes
 369   if (GetDebug()>2)
 370     cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
 371        << ogime->TreeD()->GetName()<<endl;
 372
 373    ogime->WriteDigits("OVERWRITE");
 374
 375   delete [] fHitMap;
 376   delete fTDList;
 377
 378   if (fHits)    fHits->Delete();
 379   if (fPadHits) fPadHits->Delete();
 380 }
 381
 382
 383 //------------------------------------------------------------------------
 384 void AliMUONDigitizer::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 385 {
 386   //
 387   // Sort the list of tracks contributing to a given digit
 388   // Only the 3 most significant tracks are acctually sorted
 389   //
 390
 391   //
 392   //  Loop over signals, only 3 times
 393   //
 394
 395   Int_t qmax;
 396   Int_t jmax;
 397   Int_t idx[3] = {-2,-2,-2};
 398   Int_t jch[3] = {-2,-2,-2};
 399   Int_t jtr[3] = {-2,-2,-2};
 400   Int_t i,j,imax;
 401
 402   if (ntr<3) imax=ntr;
 403   else imax=3;
 404   for(i=0;i<imax;i++){
 405     qmax=0;
 406     jmax=0;
 407
 408     for(j=0;j<ntr;j++){
 409
 410       if((i == 1 && j == idx[i-1])
 411         ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 412
 413       if(charges[j] > qmax) {
 414        qmax = charges[j];
 415        jmax=j;
 416       }
 417     }
 418
 419     if(qmax > 0) {
 420       idx[i]=jmax;
 421       jch[i]=charges[jmax];
 422       jtr[i]=tracks[jmax];
 423     }
 424
 425   }
 426
 427   for(i=0;i<3;i++){
 428     if (jtr[i] == -2) {
 429       charges[i]=0;
 430       tracks[i]=0;
 431     } else {
 432       charges[i]=jch[i];
 433       tracks[i]=jtr[i];
 434     }
 435   }
 436 }