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