MUON/AliMUONDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.3  2001/10/04 20:01:54  jchudoba
  19 changes for TTask implementation, some other small editing
  20
  21 Revision 1.2  2001/07/28 10:46:04  hristov
  22 AliRunDigitizer.h included; typos corrected
  23
  24 Revision 1.1  2001/07/27 15:41:01  jchudoba
  25 merging/digitization classes
  26
  27 */
  28 #include <TTree.h>
  29 #include <TVector.h>
  30 #include <TObjArray.h>
  31 #include <TFile.h>
  32 #include <TDirectory.h>
  33 #include <iostream.h>
  34
  35 #include "AliMUONDigitizer.h"
  36 #include "AliMUONConstants.h"
  37 #include "AliMUONChamber.h"
  38 #include "AliHitMap.h"
  39 #include "AliMUONHitMapA1.h"
  40 #include "AliMUON.h"
  41 #include "AliMUONHit.h"
  42 #include "AliMUONPadHit.h"
  43 #include "AliMUONDigit.h"
  44 #include "AliMUONTransientDigit.h"
  45 #include "AliRun.h"
  46 #include "AliPDG.h"
  47 #include "AliRunDigitizer.h"
  48
  49 ClassImp(AliMUONDigitizer)
  50
  51 //___________________________________________
  52 AliMUONDigitizer::AliMUONDigitizer() :AliDigitizer()
  53 {
  54 // Default ctor - don't use it
  55   ;
  56 }
  57
  58 //___________________________________________
  59 AliMUONDigitizer::AliMUONDigitizer(AliRunDigitizer* manager)
  60     :AliDigitizer(manager)
  61 {
  62 // ctor which should be used
  63   fHitMap  = 0;
  64   fTDList  = 0;
  65   fHits    = 0;
  66   fPadHits = 0;
  67   fTrList  = 0;
  68   fAddress = 0;
  69   fDebug   = 0;
  70   if (GetDebug()>2)
  71     cerr<<"AliMUONDigitizer::AliMUONDigitizer"
  72         <<"(AliRunDigitizer* manager) was processed"<<endl;
  73 }
  74
  75 //------------------------------------------------------------------------
  76 AliMUONDigitizer::~AliMUONDigitizer()
  77 {
  78 // Destructor
  79   if (fHits)       delete fHits;
  80   if (fPadHits)    delete fPadHits;
  81   if (fHitMap)     delete fHitMap;
  82   if (fTDList)     delete fTDList;
  83   if (fTrList)     delete fTrList;
  84   if (fAddress)    delete fAddress;
  85 }
  86
  87 //------------------------------------------------------------------------
  88 Bool_t AliMUONDigitizer::Exists(const AliMUONPadHit *mergable)
  89 {
  90   AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
  91   return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
  92 }
  93
  94 //------------------------------------------------------------------------
  95 void AliMUONDigitizer::Update(AliMUONPadHit *mergable)
  96 {
  97   AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
  98   AliMUONTransientDigit* pdigit;
  99   Int_t ipx      = padhit->PadX();        // pad number on X
 100   Int_t ipy      = padhit->PadY();        // pad number on Y
 101   Int_t iqpad    = Int_t(padhit->QPad()); // charge per pad
 102
 103   pdigit = (AliMUONTransientDigit*) fHitMap[fNch]->GetHit(ipx, ipy);
 104   // update charge
 105   //
 106   (*pdigit).AddSignal(iqpad);
 107   (*pdigit).AddPhysicsSignal(iqpad);
 108   // update list of tracks
 109   //
 110   TObjArray* fTrList = (TObjArray*)pdigit->TrackList();
 111   Int_t lastEntry = fTrList->GetLast();
 112   TVector *pTrack = (TVector*)fTrList->At(lastEntry);
 113   TVector &ptrk   = *pTrack;
 114   TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
 115   Int_t lastTrack = Int_t(ptrk(0));
 116
 117   if (trinfo(0) != kBgTag) {
 118     if (lastTrack == fTrack) {
 119       Int_t lastCharge = Int_t(ptrk(1));
 120       lastCharge += iqpad;
 121       fTrList->RemoveAt(lastEntry);
 122       trinfo(1) = lastCharge;
 123       fTrList->AddAt(&trinfo,lastEntry);
 124     } else {
 125       fTrList->Add(&trinfo);
 126     }
 127   } else {
 128     if (lastTrack != -1) fTrList->Add(&trinfo);
 129   }
 130 }
 131
 132 //------------------------------------------------------------------------
 133 void AliMUONDigitizer::CreateNew(AliMUONPadHit *mergable)
 134 {
 135 // Create new AliMUONTransientDigit and add it to the fTDList
 136
 137   AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
 138   AliMUONTransientDigit* pdigit;
 139
 140   Int_t ipx      = padhit->PadX();       // pad number on X
 141   Int_t ipy      = padhit->PadY();       // pad number on Y
 142   fTDList->AddAtAndExpand(
 143     new AliMUONTransientDigit(fNch,fDigits),fCounter);
 144   fHitMap[fNch]->SetHit(ipx, ipy, fCounter);
 145   fCounter++;
 146   pdigit = (AliMUONTransientDigit*)fTDList->At(fTDList->GetLast());
 147   // list of tracks
 148   TObjArray *fTrList = (TObjArray*)pdigit->TrackList();
 149   TVector &trinfo    = *((TVector*) (*fAddress)[fCountadr-1]);
 150   fTrList->Add(&trinfo);
 151 }
 152
 153
 154 //------------------------------------------------------------------------
 155 Bool_t AliMUONDigitizer::Init()
 156 {
 157 // Initialization
 158
 159   fHits     = new TClonesArray("AliMUONHit",1000);
 160   fPadHits  = new TClonesArray("AliMUONPadHit",1000);
 161   return kTRUE;
 162 }
 163
 164
 165 //------------------------------------------------------------------------
 166 //void AliMUONDigitizer::Digitize()
 167 void AliMUONDigitizer::Exec(Option_t* option)
 168 {
 169
 170   TString optionString = option;
 171   if (optionString.Data() == "deb") {
 172     cout<<"AliMUONDigitizer::Exec: called with option deb "<<endl;
 173     fDebug = 3;
 174   }
 175   AliMUONChamber*   iChamber;
 176   AliSegmentation*  segmentation;
 177
 178   if (GetDebug()>2) cerr<<"   AliMUONDigitizer::Digitize() starts"<<endl;
 179   fTDList = new TObjArray;
 180   if(!fAddress) fAddress = new TClonesArray("TVector",10000);
 181
 182   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 183   if (!pMUON) {
 184     cerr<<"AliMUONDigitizer::Digitize Error:"
 185         <<" module MUON not found in the input file"<<endl;
 186     return;
 187   }
 188   pMUON->MakeBranchInTreeD(fManager->GetTreeD());
 189   fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
 190   for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {fHitMap[i] = 0;}
 191
 192   //
 193   // loop over cathodes
 194   //
 195
 196   fCountadr = 0;
 197   for (int icat = 0; icat < 2; icat++) {
 198     fCounter = 0;
 199     Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
 200     for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
 201       iChamber = &(pMUON->Chamber(i));
 202 //      if (!(iChamber->Nsec() == 1 && icat == 1)) {
 203         segmentation = iChamber->SegmentationModel(icat+1);
 204         fHitMap[i] = new AliMUONHitMapA1(segmentation, fTDList);
 205         nmuon[i] = 0;
 206 //      }
 207     }
 208
 209
 210 // Loop over files to digitize
 211
 212     for (Int_t inputFile=0; inputFile<fManager->GetNinputs();
 213          inputFile++) {
 214 // Connect MUON branches
 215
 216       TBranch *branchHits = 0;
 217       TBranch *branchPadHits = 0;
 218       TTree *treeH = fManager->GetInputTreeH(inputFile);
 219       if (GetDebug()>2) {
 220         cerr<<"   inputFile , cathode = "<<inputFile<<" "
 221             <<icat<<endl;
 222         cerr<<"   treeH, fHits "<<treeH<<" "<<fHits<<endl;
 223       }
 224       if (treeH && fHits) {
 225         branchHits = treeH->GetBranch("MUON");
 226         if (branchHits) {
 227           fHits->Clear();
 228           branchHits->SetAddress(&fHits);
 229         }
 230         else
 231           Error("Digitize","branch MUON was not found");
 232       }
 233       if (GetDebug()>2) cerr<<"   branchHits = "<<branchHits<<endl;
 234
 235       if (treeH && fPadHits) {
 236         branchPadHits = treeH->GetBranch("MUONCluster");
 237         if (branchPadHits)
 238           branchPadHits->SetAddress(&fPadHits);
 239         else
 240           Error("Digitize","branch MUONCluster was not found");
 241       }
 242       if (GetDebug()>2) cerr<<"   branchPadHits = "<<branchPadHits<<endl;
 243
 244 //
 245 //   Loop over tracks
 246 //
 247
 248       Int_t ntracks = (Int_t) treeH->GetEntries();
 249
 250       for (fTrack = 0; fTrack < ntracks; fTrack++) {
 251         if (GetDebug()>2) cerr<<"   fTrack = "<<fTrack<<endl;
 252         fHits->Clear();
 253         fPadHits->Clear();
 254         branchHits->GetEntry(fTrack);
 255         branchPadHits->GetEntry(fTrack);
 256
 257 //
 258 //   Loop over hits
 259
 260         AliMUONHit* mHit;
 261         for(int i = 0; i < fHits->GetEntriesFast(); ++i)
 262         {
 263           mHit   = (AliMUONHit*) (*fHits)[i];
 264           fNch = mHit->Chamber()-1;  // chamber number
 265           iChamber = &(pMUON->Chamber(fNch));
 266           if (fNch > AliMUONConstants::NCh()-1) {
 267             cerr<<"AliMUONDigitizer: ERROR: "
 268                 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
 269                 <<fNch<<", "<< AliMUONConstants::NCh()<<endl;
 270             return;
 271           }
 272 //
 273 // Loop over pad hits
 274           for (AliMUONPadHit* mPad =
 275                  (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
 276                mPad;
 277                mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHits))
 278           {
 279             Int_t cathode  = mPad->Cathode();      // cathode number
 280             Int_t ipx      = mPad->PadX();         // pad number on X
 281             Int_t ipy      = mPad->PadY();         // pad number on Y
 282             Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
 283             if (cathode != (icat+1)) continue;
 284
 285             new((*fAddress)[fCountadr++]) TVector(2);
 286
 287             TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
 288             Int_t mask = fManager->GetMask(inputFile);
 289             trinfo(0) = (Float_t)(fTrack + mask);  // tag background
 290 //                  trinfo(0) = (Float_t)fTrack;
 291             if (inputFile == 0) {
 292               trinfo(1) = (Float_t)iqpad;
 293             } else {
 294               trinfo(1) = kBgTag;
 295             }
 296             fDigits[0] = ipx;
 297             fDigits[1] = ipy;
 298             fDigits[2] = icat;
 299             fDigits[3] = iqpad;
 300             if (inputFile == 0) {
 301               fDigits[4] = iqpad;
 302             } else {
 303               fDigits[4] = 0;
 304             }
 305             if (mHit->Particle() == kMuonPlus ||
 306                 mHit->Particle() == kMuonMinus) {
 307               fDigits[5] = (mPad->HitNumber()) + mask;
 308             } else fDigits[5] = -1;
 309
 310             // build the list of fired pads and update the info,
 311             // fDigits is input for Update(mPad)
 312
 313             if (!Exists(mPad)) {
 314               CreateNew(mPad);
 315             } else {
 316               Update(mPad);
 317             } //  end if Exists(mPad)
 318           } //end loop over clusters
 319         } // hit loop
 320       } // track loop
 321     } // end file loop
 322     if (GetDebug()>2) cerr<<"END OF FILE LOOP"<<endl;
 323
 324     Int_t tracks[MAXTRACKS];
 325     Int_t charges[MAXTRACKS];
 326     Int_t nentries = fTDList->GetEntriesFast();
 327
 328     for (Int_t nent = 0; nent < nentries; nent++) {
 329       AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
 330       if (address == 0) continue;
 331       Int_t ich = address->Chamber();
 332       Int_t   q = address->Signal();
 333       iChamber = &(pMUON->Chamber(ich));
 334 //
 335 //  Digit Response (noise, threshold, saturation, ...)
 336       AliMUONResponse * response = iChamber->ResponseModel();
 337       q = response->DigitResponse(q);
 338
 339       if (!q) continue;
 340
 341       fDigits[0] = address->PadX();
 342       fDigits[1] = address->PadY();
 343       fDigits[2] = address->Cathode();
 344       fDigits[3] = q;
 345       fDigits[4] = address->Physics();
 346       fDigits[5] = address->Hit();
 347
 348       TObjArray* fTrList = (TObjArray*)address->TrackList();
 349       Int_t nptracks = fTrList->GetEntriesFast();
 350
 351       // this was changed to accomodate the real number of tracks
 352
 353       if (nptracks > MAXTRACKS) {
 354         if (GetDebug() >0) {
 355           cerr<<"AliMUONDigitizer: nptracks > 10 "<<nptracks;
 356           cerr<<"reset to max value "<<MAXTRACKS<<endl;
 357         }
 358         nptracks = MAXTRACKS;
 359       }
 360       if (nptracks > 2 && GetDebug() >2) {
 361         cerr<<"AliMUONDigitizer: nptracks > 2 "<<nptracks<<endl;
 362         printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
 363       }
 364       for (Int_t tr = 0; tr < nptracks; tr++) {
 365         TVector *ppP = (TVector*)fTrList->At(tr);
 366         if(!ppP ) {
 367           cerr<<"Error: ppP = "<<ppP<<endl;
 368           return;
 369         }
 370         TVector &pp  = *ppP;
 371         tracks[tr]   = Int_t(pp(0));
 372         charges[tr]  = Int_t(pp(1));
 373       }      //end loop over list of tracks for one pad
 374       // Sort list of tracks according to charge
 375       if (nptracks > 1) {
 376         SortTracks(tracks,charges,nptracks);
 377       }
 378       if (nptracks < MAXTRACKS ) {
 379         for (Int_t i = nptracks; i < MAXTRACKS; i++) {
 380           tracks[i]  = 0;
 381           charges[i] = 0;
 382         }
 383       }
 384
 385       // fill digits
 386       pMUON->AddDigits(ich,tracks,charges,fDigits);
 387       // delete fTrList;
 388     }
 389
 390     fManager->GetTreeD()->Fill();
 391
 392     pMUON->ResetDigits();  //
 393     fTDList->Clear();
 394
 395
 396     for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
 397       if (fHitMap[ii]) {
 398         delete fHitMap[ii];
 399         fHitMap[ii] = 0;
 400       }
 401     }
 402     delete [] nmuon;
 403   } //end loop over cathodes
 404   if (GetDebug()>2)
 405     cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
 406         <<fManager->GetTreeD()->GetName()<<endl;
 407   fManager->GetTreeD()->Write(0,TObject::kOverwrite);
 408   delete [] fHitMap;
 409   delete fTDList;
 410
 411   if (fAddress)    fAddress->Delete();
 412   if (fHits)    fHits->Delete();
 413   if (fPadHits) fPadHits->Delete();
 414 }
 415
 416
 417 //------------------------------------------------------------------------
 418 void AliMUONDigitizer::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 419 {
 420   //
 421   // Sort the list of tracks contributing to a given digit
 422   // Only the 3 most significant tracks are acctually sorted
 423   //
 424
 425   //
 426   //  Loop over signals, only 3 times
 427   //
 428
 429   Int_t qmax;
 430   Int_t jmax;
 431   Int_t idx[3] = {-2,-2,-2};
 432   Int_t jch[3] = {-2,-2,-2};
 433   Int_t jtr[3] = {-2,-2,-2};
 434   Int_t i,j,imax;
 435
 436   if (ntr<3) imax=ntr;
 437   else imax=3;
 438   for(i=0;i<imax;i++){
 439     qmax=0;
 440     jmax=0;
 441
 442     for(j=0;j<ntr;j++){
 443
 444       if((i == 1 && j == idx[i-1])
 445          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 446
 447       if(charges[j] > qmax) {
 448         qmax = charges[j];
 449         jmax=j;
 450       }
 451     }
 452
 453     if(qmax > 0) {
 454       idx[i]=jmax;
 455       jch[i]=charges[jmax];
 456       jtr[i]=tracks[jmax];
 457     }
 458
 459   }
 460
 461   for(i=0;i<3;i++){
 462     if (jtr[i] == -2) {
 463       charges[i]=0;
 464       tracks[i]=0;
 465     } else {
 466       charges[i]=jch[i];
 467       tracks[i]=jtr[i];
 468     }
 469   }
 470 }