MUON/AliMUONDigitizerv1.cxx

   1
   2
   3 #include <Riostream.h>
   4 #include <TDirectory.h>
   5 #include <TFile.h>
   6 #include <TObjArray.h>
   7 #include <TPDGCode.h>
   8 #include <TTree.h>
   9
  10 #include "AliMUON.h"
  11 #include "AliMUONChamber.h"
  12 #include "AliMUONConstants.h"
  13 #include "AliMUONDigit.h"
  14 #include "AliMUONDigitizerv1.h"
  15 #include "AliMUONHit.h"
  16 #include "AliMUONHitMapA1.h"
  17 #include "AliMUONPadHit.h"
  18 #include "AliMUONTransientDigit.h"
  19 #include "AliRun.h"
  20 #include "AliRunDigitizer.h"
  21
  22 ClassImp(AliMUONDigitizerv1)
  23
  24 //___________________________________________
  25 AliMUONDigitizerv1::AliMUONDigitizerv1() :AliDigitizer()
  26 {
  27 // Default ctor - don't use it
  28   fHitMap = 0;
  29   fTDList = 0;
  30   if (GetDebug()>2)
  31     cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
  32         <<"(AliRunDigitizer* manager) was processed"<<endl;
  33 }
  34
  35 //___________________________________________
  36 AliMUONDigitizerv1::AliMUONDigitizerv1(AliRunDigitizer* manager)
  37     :AliDigitizer(manager)
  38 {
  39 // ctor which should be used
  40   fHitMap  = 0;
  41   fTDList  = 0;
  42   fDebug   = 0;
  43   fHits = new TClonesArray("AliMUONHit",1000);
  44   if (GetDebug()>2)
  45     cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
  46         <<"(AliRunDigitizer* manager) was processed"<<endl;
  47 }
  48
  49 //------------------------------------------------------------------------
  50 AliMUONDigitizerv1::~AliMUONDigitizerv1()
  51 {
  52 // Destructor
  53   delete fHits;
  54 }
  55
  56 //------------------------------------------------------------------------
  57 void AliMUONDigitizerv1::AddTransientDigit(AliMUONTransientDigit * mTD)
  58 {
  59   // Choosing the maping of the cathode plane of the chamber:
  60   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  61   fTDList->AddAtAndExpand(mTD, fTDCounter);
  62   fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
  63   fTDCounter++;
  64 }
  65
  66 //------------------------------------------------------------------------
  67 Bool_t AliMUONDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD)
  68 {
  69   // Choosing the maping of the cathode plane of the chamber:
  70   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  71   return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
  72 }
  73
  74 //------------------------------------------------------------------------
  75 Bool_t AliMUONDigitizerv1::Init()
  76 {
  77 // Initialization
  78   return kTRUE;
  79 }
  80
  81 //------------------------------------------------------------------------
  82 void AliMUONDigitizerv1::UpdateTransientDigit(Int_t track, AliMUONTransientDigit * mTD)
  83 {
  84   // Choosing the maping of the cathode plane of the chamber:
  85   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  86   AliMUONTransientDigit *pdigit =
  87     static_cast<AliMUONTransientDigit*>(fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()));
  88   // update charge
  89   //
  90   Int_t iqpad  = mTD->Signal();        // charge per pad
  91   pdigit->AddSignal(iqpad);
  92   pdigit->AddPhysicsSignal(iqpad);
  93   // update list of tracks
  94   //
  95   Int_t charge;
  96   track=+ fMask;
  97   if (fSignal) charge = iqpad;
  98   else         charge = kBgTag;
  99   pdigit->UpdateTrackList(track,charge);
 100 }
 101
 102
 103 //--------------------------------------------------------------------------
 104 void AliMUONDigitizerv1::MakeTransientDigit(Int_t track, Int_t iHit, AliMUONHit * mHit)
 105 {
 106   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 107   if (!pMUON) {
 108     cerr<<"AliMUONDigitizerv1::Digitize Error:"
 109         <<" module MUON not found in the input file"<<endl;
 110   }
 111   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit starts"<<endl;
 112   Int_t   ichamber = mHit->Chamber()-1;
 113   AliMUONChamber & chamber = pMUON->Chamber(ichamber);
 114   Float_t xhit = mHit->X();
 115   Float_t yhit = mHit->Y();
 116   Float_t zhit = mHit->Z();
 117   Float_t eloss= mHit->Eloss();
 118   Float_t tof  = mHit->Age();
 119   // Variables for chamber response from AliMUONChamber::DisIntegration
 120   Float_t newdigit[6][500];  // Pad information
 121   Int_t nnew=0;              // Number of touched Pads per hit
 122   Int_t digits[6];
 123
 124   //
 125   // Calls the charge disintegration method of the current chamber
 126   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit calling AliMUONChamber::DisIngtegration starts"<<endl;
 127   chamber.DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newdigit);
 128   // Creating a new TransientDigits from hit
 129   for(Int_t iTD=0; iTD<nnew; iTD++) {
 130     digits[0] = Int_t(newdigit[1][iTD]);  // Padx of the Digit
 131     digits[1] = Int_t(newdigit[2][iTD]);  // Pady of the Digit
 132     digits[2] = Int_t(newdigit[5][iTD]);  // Cathode plane
 133     digits[3] = Int_t(newdigit[3][iTD]);  // Induced charge in the Pad
 134     if (fSignal) digits[4] = Int_t(newdigit[3][iTD]);
 135     else         digits[4] = 0;
 136     digits[5] = iHit+fMask;    // Hit number in the list
 137     if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit " <<
 138                         "PadX "<< digits[0] << " " <<
 139                         "PadY "<< digits[1] << " " <<
 140                         "Plane " << digits[2] << " " <<
 141                         "Charge " << digits[3] <<" " <<
 142                         "Hit " << digits[5] << endl;
 143     // list of tracks
 144     Int_t charge;
 145     track += fMask;
 146     if (fSignal) charge = digits[3];
 147     else         charge = kBgTag;
 148     if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
 149     AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
 150     mTD->AddToTrackList(track,charge);
 151
 152     if (!ExistTransientDigit(mTD)) {
 153       AddTransientDigit(mTD);
 154       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Adding TransientDigit"<<endl;
 155     }
 156     else {
 157       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit updating TransientDigit"<<endl;
 158       UpdateTransientDigit(track, mTD);
 159       delete mTD;
 160     }
 161   }
 162 }
 163 //-----------------------------------------------------------------------
 164 void AliMUONDigitizerv1::Exec(Option_t* option)
 165 {
 166   TString optionString = option;
 167   if (optionString.Data() == "deb") {
 168     cout<<"AliMUONDigitizerv1::Exec: called with option deb "<<endl;
 169     fDebug = 3;
 170   }
 171
 172   AliMUONChamber*   chamber;
 173   AliSegmentation*  c1Segmentation; //Cathode plane c1 of the chamber
 174   AliSegmentation*  c2Segmentation; //Cathode place c2 of the chamber
 175
 176   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() starts"<<endl;
 177   fTDList = new TObjArray;
 178
 179   // Getting Module MUON
 180   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 181   if (!pMUON) {
 182     cerr<<"AliMUONDigitizerv1::Digitize Error:"
 183         <<" module MUON not found in the input file"<<endl;
 184     return;
 185   }
 186   // New branch for MUON digit in the tree of digits
 187   pMUON->MakeBranchInTreeD(fManager->GetTreeD());
 188
 189   // Array of pointer of the AliMUONHitMapA1:
 190   //  two HitMaps per chamber, or one HitMap per cahtode plane
 191   fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
 192
 193   //Loop over chambers for the definition AliMUONHitMap
 194   for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
 195     chamber = &(pMUON->Chamber(i));
 196     c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
 197     fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
 198     c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
 199     fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
 200   }
 201
 202 // Loop over files to merge and to digitize
 203     fSignal = kTRUE;
 204     for (Int_t inputFile=0; inputFile<fManager->GetNinputs(); inputFile++) {
 205       // Connect MUON Hit branch
 206       if (inputFile > 0 ) fSignal = kFALSE;
 207       TBranch *branchHits = 0;
 208       TTree *treeH = fManager->GetInputTreeH(inputFile);
 209       if (GetDebug()>2) {
 210         cerr<<"AliMUONDigitizerv1::Exec inputFile is "<<inputFile<<" "<<endl;
 211         cerr<<"AliMUONDigitizerv1::Exec treeH, fHits "<<treeH<<" "<<fHits<<endl;
 212       }
 213       if (treeH && fHits) {
 214         branchHits = treeH->GetBranch("MUON");
 215         if (branchHits) {
 216           fHits->Clear();
 217           branchHits->SetAddress(&fHits);
 218         }
 219         else
 220           Error("Exec","branch MUON was not found");
 221       }
 222       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec branchHits = "<<branchHits<<endl;
 223
 224       fMask = fManager->GetMask(inputFile);
 225       //
 226       // Loop over tracks
 227       Int_t itrack;
 228       Int_t ntracks = (Int_t) treeH->GetEntries();
 229       for (itrack = 0; itrack < ntracks; itrack++) {
 230         if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec itrack = "<<itrack<<endl;
 231         fHits->Clear();
 232         branchHits->GetEntry(itrack);
 233         //
 234         //  Loop over hits
 235         Int_t ihit, ichamber;
 236         AliMUONHit* mHit;
 237         for(ihit = 0; ihit < fHits->GetEntriesFast(); ihit++) {
 238           mHit = static_cast<AliMUONHit*>(fHits->At(ihit));
 239           ichamber = mHit->Chamber()-1;  // chamber number
 240           if (ichamber > AliMUONConstants::NCh()-1) {
 241             cerr<<"AliMUONDigitizer: ERROR: "
 242                 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
 243                 <<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
 244             return;
 245           }
 246           chamber = &(pMUON->Chamber(ichamber));
 247           //
 248           //Dumping Hit content:
 249           if (GetDebug()>2) {
 250             cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
 251               ihit << " " <<
 252               mHit->Chamber() << " " <<
 253               mHit->X() << " " <<
 254               mHit->Y() << " " <<
 255               mHit->Z() << " " <<
 256               mHit->Eloss() << " " << endl;
 257           }
 258           //
 259           // Inititializing Correlation
 260           chamber->ChargeCorrelationInit();
 261           if (ichamber < AliMUONConstants::NTrackingCh()) {
 262             // Tracking Chamber
 263             // Initialize hit position (cursor) in the segmentation model
 264             chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
 265           } else {
 266             // Trigger Chamber
 267           }
 268           MakeTransientDigit(itrack, ihit, mHit);
 269         } // hit loop
 270       } // track loop
 271     } // end file loop
 272     if (GetDebug()>2) cerr<<"AliMUONDigitizer::Exec End of hits, track and file loops"<<endl;
 273
 274
 275     // Loop on cathodes
 276     Int_t icat;
 277     for(icat=0; icat<2; icat++) {
 278       //
 279       // Filling Digit List
 280       Int_t tracks[kMAXTRACKS];
 281       Int_t charges[kMAXTRACKS];
 282       Int_t nentries = fTDList->GetEntriesFast();
 283       Int_t digits[6];
 284       for (Int_t nent = 0; nent < nentries; nent++) {
 285         AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
 286         if (address == 0) continue;
 287         Int_t ich = address->Chamber();
 288         Int_t   q = address->Signal();
 289         chamber = &(pMUON->Chamber(ich));
 290         //
 291         //  Digit Response (noise, threshold, saturation, ...)
 292         AliMUONResponse * response = chamber->ResponseModel();
 293         q = response->DigitResponse(q,address);
 294
 295         if (!q) continue;
 296
 297         digits[0] = address->PadX();
 298         digits[1] = address->PadY();
 299         digits[2] = address->Cathode()-1;
 300         digits[3] = q;
 301         digits[4] = address->Physics();
 302         digits[5] = address->Hit();
 303
 304         Int_t nptracks = address->GetNTracks();
 305
 306         if (nptracks > kMAXTRACKS) {
 307           if (GetDebug() >0) {
 308             cerr<<"AliMUONDigitizer:Exec  nptracks > 10 "<<nptracks;
 309             cerr<<"reset to max value "<<kMAXTRACKS<<endl;
 310           }
 311           nptracks = kMAXTRACKS;
 312         }
 313         if (nptracks > 2 && GetDebug() >2) {
 314           cerr<<"AliMUONDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
 315           //    printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
 316         }
 317         for (Int_t tr = 0; tr < nptracks; tr++) {
 318           tracks[tr]   = address->GetTrack(tr);
 319           charges[tr]  = address->GetCharge(tr);
 320         }      //end loop over list of tracks for one pad
 321         // Sort list of tracks according to charge
 322         if (nptracks > 1) {
 323           SortTracks(tracks,charges,nptracks);
 324         }
 325         if (nptracks < kMAXTRACKS ) {
 326           for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
 327             tracks[i]  = 0;
 328             charges[i] = 0;
 329           }
 330         }
 331
 332         // fill digits
 333         if (GetDebug()>2) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
 334         if ( digits[2] == icat ) pMUON->AddDigits(ich,tracks,charges,digits);
 335       }
 336       fManager->GetTreeD()->Fill();
 337       pMUON->ResetDigits();  //
 338     } // end loop cathode
 339     fTDList->Delete();
 340
 341     for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
 342       if (fHitMap[ii]) {
 343         delete fHitMap[ii];
 344         fHitMap[ii] = 0;
 345       }
 346     }
 347
 348     if (GetDebug()>2)
 349       cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
 350           <<fManager->GetTreeD()->GetName()<<endl;
 351     fManager->GetTreeD()->Write(0,TObject::kOverwrite);
 352     delete [] fHitMap;
 353     delete fTDList;
 354
 355     if (fHits)    fHits->Clear();
 356 }
 357
 358 //------------------------------------------------------------------------
 359 void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 360 {
 361   //
 362   // Sort the list of tracks contributing to a given digit
 363   // Only the 3 most significant tracks are acctually sorted
 364   //
 365
 366   //
 367   //  Loop over signals, only 3 times
 368   //
 369
 370   Int_t qmax;
 371   Int_t jmax;
 372   Int_t idx[3] = {-2,-2,-2};
 373   Int_t jch[3] = {-2,-2,-2};
 374   Int_t jtr[3] = {-2,-2,-2};
 375   Int_t i,j,imax;
 376
 377   if (ntr<3) imax=ntr;
 378   else imax=3;
 379   for(i=0;i<imax;i++){
 380     qmax=0;
 381     jmax=0;
 382
 383     for(j=0;j<ntr;j++){
 384
 385       if((i == 1 && j == idx[i-1])
 386          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 387
 388       if(charges[j] > qmax) {
 389         qmax = charges[j];
 390         jmax=j;
 391       }
 392     }
 393
 394     if(qmax > 0) {
 395       idx[i]=jmax;
 396       jch[i]=charges[jmax];
 397       jtr[i]=tracks[jmax];
 398     }
 399
 400   }
 401
 402   for(i=0;i<3;i++){
 403     if (jtr[i] == -2) {
 404       charges[i]=0;
 405       tracks[i]=0;
 406     } else {
 407       charges[i]=jch[i];
 408       tracks[i]=jtr[i];
 409     }
 410   }
 411 }