MUON/AliMUONDigitizerv2.cxx

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