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