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