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