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