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