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