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->Delete();
 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     // Filling Digit List
 276     Int_t tracks[kMAXTRACKS];
 277     Int_t charges[kMAXTRACKS];
 278     Int_t nentries = fTDList->GetEntriesFast();
 279     Int_t digits[6];
 280     for (Int_t nent = 0; nent < nentries; nent++) {
 281       AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
 282       if (address == 0) continue;
 283       Int_t ich = address->Chamber();
 284       Int_t   q = address->Signal();
 285       chamber = &(pMUON->Chamber(ich));
 286       //
 287       //  Digit Response (noise, threshold, saturation, ...)
 288       AliMUONResponse * response = chamber->ResponseModel();
 289       q = response->DigitResponse(q,address);
 290
 291       if (!q) continue;
 292
 293       digits[0] = address->PadX();
 294       digits[1] = address->PadY();
 295       digits[2] = address->Cathode();
 296       digits[3] = q;
 297       digits[4] = address->Physics();
 298       digits[5] = address->Hit();
 299
 300       Int_t nptracks = address->GetNTracks();
 301
 302       if (nptracks > kMAXTRACKS) {
 303         if (GetDebug() >0) {
 304           cerr<<"AliMUONDigitizer:Exec  nptracks > 10 "<<nptracks;
 305           cerr<<"reset to max value "<<kMAXTRACKS<<endl;
 306         }
 307         nptracks = kMAXTRACKS;
 308       }
 309       if (nptracks > 2 && GetDebug() >2) {
 310         cerr<<"AliMUONDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
 311         //      printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
 312       }
 313       for (Int_t tr = 0; tr < nptracks; tr++) {
 314         tracks[tr]   = address->GetTrack(tr);
 315         charges[tr]  = address->GetCharge(tr);
 316       }      //end loop over list of tracks for one pad
 317       // Sort list of tracks according to charge
 318       if (nptracks > 1) {
 319         SortTracks(tracks,charges,nptracks);
 320       }
 321       if (nptracks < kMAXTRACKS ) {
 322         for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
 323           tracks[i]  = 0;
 324           charges[i] = 0;
 325         }
 326       }
 327
 328       // fill digits
 329       if (GetDebug()>2) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
 330       pMUON->AddDigits(ich,tracks,charges,digits);
 331     }
 332     fManager->GetTreeD()->Fill();
 333     pMUON->ResetDigits();  //
 334     fTDList->Delete();
 335
 336     for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
 337       if (fHitMap[ii]) {
 338         delete fHitMap[ii];
 339         fHitMap[ii] = 0;
 340       }
 341     }
 342
 343     if (GetDebug()>2)
 344       cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
 345           <<fManager->GetTreeD()->GetName()<<endl;
 346     fManager->GetTreeD()->Write(0,TObject::kOverwrite);
 347     delete [] fHitMap;
 348     delete fTDList;
 349
 350     if (fHits)    fHits->Clear();
 351 }
 352
 353 //------------------------------------------------------------------------
 354 void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 355 {
 356   //
 357   // Sort the list of tracks contributing to a given digit
 358   // Only the 3 most significant tracks are acctually sorted
 359   //
 360
 361   //
 362   //  Loop over signals, only 3 times
 363   //
 364
 365   Int_t qmax;
 366   Int_t jmax;
 367   Int_t idx[3] = {-2,-2,-2};
 368   Int_t jch[3] = {-2,-2,-2};
 369   Int_t jtr[3] = {-2,-2,-2};
 370   Int_t i,j,imax;
 371
 372   if (ntr<3) imax=ntr;
 373   else imax=3;
 374   for(i=0;i<imax;i++){
 375     qmax=0;
 376     jmax=0;
 377
 378     for(j=0;j<ntr;j++){
 379
 380       if((i == 1 && j == idx[i-1])
 381          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 382
 383       if(charges[j] > qmax) {
 384         qmax = charges[j];
 385         jmax=j;
 386       }
 387     }
 388
 389     if(qmax > 0) {
 390       idx[i]=jmax;
 391       jch[i]=charges[jmax];
 392       jtr[i]=tracks[jmax];
 393     }
 394
 395   }
 396
 397   for(i=0;i<3;i++){
 398     if (jtr[i] == -2) {
 399       charges[i]=0;
 400       tracks[i]=0;
 401     } else {
 402       charges[i]=jch[i];
 403       tracks[i]=jtr[i];
 404     }
 405   }
 406 }