MUON/AliMUONMerger.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.10  2003/01/14 10:50:19  alibrary
  19 Cleanup of STEER coding conventions
  20
  21 Revision 1.9  2002/03/13 07:04:11  jchudoba
  22 Connect only MUON branches when reading the event to speed up digitisation.
  23
  24 Revision 1.8  2002/02/22 12:14:21  morsch
  25 Validate pad hit before digitization.
  26
  27 Revision 1.7  2001/11/22 11:15:41  jchudoba
  28 Proper deletion of arrays (thanks to Rene Brun)
  29
  30 Revision 1.6  2001/11/02 12:55:45  jchudoba
  31 cleanup of the code, add const to Get methods
  32
  33 Revision 1.5  2001/10/31 16:35:07  jchudoba
  34 some functionality move to AliMUONTransientDigit class
  35
  36 Revision 1.4  2001/03/20 13:36:11  egangler
  37 TFile memory leak and "too many files open" problem solved (?):
  38 the backgroud file is open only once forever, and not once for each event.
  39
  40 Revision 1.3  2001/03/05 23:57:44  morsch
  41 Writing of digit tree moved to macro.
  42
  43 Revision 1.2  2001/03/05 08:40:25  morsch
  44 Method SortTracks(..) imported from AliMUON.
  45
  46 Revision 1.1  2001/02/02 14:11:53  morsch
  47 AliMUONMerger prototype to be called by the merge manager.
  48
  49 */
  50
  51 #include <TDirectory.h>
  52 #include <TFile.h>
  53 #include <TObjArray.h>
  54 #include <TPDGCode.h>
  55 #include <TTree.h>
  56
  57 #include "AliHitMap.h"
  58 #include "AliMUON.h"
  59 #include "AliMUONChamber.h"
  60 #include "AliMUONConstants.h"
  61 #include "AliMUONDigit.h"
  62 #include "AliMUONHit.h"
  63 #include "AliMUONHitMapA1.h"
  64 #include "AliMUONMerger.h"
  65 #include "AliMUONPadHit.h"
  66 #include "AliMUONTransientDigit.h"
  67 #include "AliRun.h"
  68
  69 ClassImp(AliMUONMerger)
  70
  71 //___________________________________________
  72 AliMUONMerger::AliMUONMerger()
  73 {
  74 // Default constructor
  75     fEvNrSig = 0;
  76     fEvNrBgr = 0;
  77     fMerge   = kDigitize;
  78     fFnBgr   = 0;
  79     fHitMap  = 0;
  80     fList    = 0;
  81     fTrH1    = 0;
  82     fHitsBgr = 0;
  83     fPadHitsBgr = 0;
  84     fHitMap     = 0;
  85     fList       = 0;
  86     fBgrFile    = 0;
  87     fDebug      = 0;
  88 }
  89
  90 //------------------------------------------------------------------------
  91 AliMUONMerger::~AliMUONMerger()
  92 {
  93 // Destructor
  94     if (fTrH1)       delete fTrH1;
  95     if (fHitsBgr)    delete fHitsBgr;
  96     if (fPadHitsBgr) delete fPadHitsBgr;
  97     if (fHitMap)     delete [] fHitMap;
  98     if (fList)       delete fList;
  99     if (fBgrFile)    delete fBgrFile;
 100 }
 101
 102 //------------------------------------------------------------------------
 103 Bool_t AliMUONMerger::Exists(const AliMUONPadHit *padhit) const
 104 {
 105 // test if the given padhit was already fired
 106     return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
 107 }
 108
 109 //------------------------------------------------------------------------
 110 void AliMUONMerger::Update(AliMUONPadHit *padhit)
 111 {
 112 // add new contribution to the fired padhit
 113     AliMUONTransientDigit *pdigit =
 114       static_cast<AliMUONTransientDigit*>(
 115       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
 116
 117     // update charge
 118     //
 119     Int_t iqpad    = padhit->QPad();        // charge per pad
 120     pdigit->AddSignal(iqpad);
 121     pdigit->AddPhysicsSignal(iqpad);
 122
 123     // update list of tracks
 124     //
 125     Int_t track, charge;
 126     if (fSignal) {
 127       track = fTrack;
 128       charge = iqpad;
 129     } else {
 130       track = kBgTag;
 131       charge = kBgTag;
 132     }
 133     pdigit->UpdateTrackList(track,charge);
 134 }
 135
 136 //------------------------------------------------------------------------
 137 void AliMUONMerger::CreateNew(AliMUONPadHit *padhit)
 138 {
 139 // add new transient digit to the list, update hit map
 140     fList->AddAtAndExpand(
 141         new AliMUONTransientDigit(fNch,fDigits),fCounter);
 142     fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
 143     AliMUONTransientDigit* pdigit =
 144       static_cast<AliMUONTransientDigit*>(fList->Last());
 145     // list of tracks
 146     Int_t track, charge;
 147     if (fSignal) {
 148       track = fTrack;
 149       charge = padhit->QPad();
 150     } else {
 151       track = kBgTag;
 152       charge = kBgTag;
 153     }
 154     pdigit->AddToTrackList(track,charge);
 155     fCounter++;
 156 }
 157
 158
 159 //------------------------------------------------------------------------
 160 void AliMUONMerger::Init()
 161 {
 162 // Initialisation
 163     // open only once the background file !!
 164     if (fMerge && !fBgrFile) fBgrFile = InitBgr();
 165 }
 166
 167
 168
 169 //------------------------------------------------------------------------
 170 TFile* AliMUONMerger::InitBgr()
 171 {
 172 // Initialise background event
 173     TFile *file = new TFile(fFnBgr);
 174 // add error checking later
 175     printf("\n AliMUONMerger has opened %s file with background event \n", fFnBgr);
 176     fHitsBgr     = new TClonesArray("AliMUONHit",1000);
 177     fPadHitsBgr  = new TClonesArray("AliMUONPadHit",1000);
 178     return file;
 179 }
 180
 181 //------------------------------------------------------------------------
 182 void AliMUONMerger::Digitise()
 183 {
 184
 185     // keep galice.root for signal and name differently the file for
 186     // background when add! otherwise the track info for signal will be lost !
 187
 188     AliMUONChamber*   iChamber;
 189     AliSegmentation*  segmentation;
 190
 191     fList = new TObjArray;
 192
 193     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 194     fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
 195     if (fMerge ) {
 196         fBgrFile->cd();
 197         //
 198         // Get Hits Tree header from file
 199         //if(fHitsBgr) fHitsBgr->Clear();     // Useless because line 327
 200         //if(fPadHitsBgr) fPadHitsBgr->Clear(); // Useless because line 328
 201         if(fTrH1) delete fTrH1;
 202         fTrH1 = 0;
 203
 204         char treeName[20];
 205         sprintf(treeName,"TreeH%d",fEvNrBgr);
 206         fTrH1 = (TTree*)gDirectory->Get(treeName);
 207         if (!fTrH1) {
 208             printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
 209         }
 210         //
 211         // Set branch addresses
 212         TBranch *branch;
 213         char branchname[20];
 214         sprintf(branchname,"%s",pMUON->GetName());
 215         if (fTrH1 && fHitsBgr) {
 216             branch = fTrH1->GetBranch(branchname);
 217             if (branch) branch->SetAddress(&fHitsBgr);
 218         }
 219         if (fTrH1 && fPadHitsBgr) {
 220             branch = fTrH1->GetBranch("MUONCluster");
 221             if (branch) branch->SetAddress(&fPadHitsBgr);
 222         }
 223     }
 224     //
 225     // loop over cathodes
 226     //
 227     fSignal = kTRUE;
 228     for (int icat = 0; icat < 2; icat++) {
 229         fCounter = 0;
 230         for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
 231             iChamber = &(pMUON->Chamber(i));
 232             if (iChamber->Nsec() == 1 && icat == 1) {
 233                 continue;
 234             } else {
 235                 segmentation = iChamber->SegmentationModel(icat+1);
 236             }
 237             fHitMap[i] = new AliMUONHitMapA1(segmentation, fList);
 238         }
 239
 240 //
 241 //   Loop over tracks
 242 //
 243
 244         TTree *treeH  = gAlice->TreeH();
 245         Int_t ntracks = (Int_t) treeH->GetEntries();
 246         treeH->SetBranchStatus("*",0); // switch off all branches
 247         treeH->SetBranchStatus("MUON*",1); // switch on only MUON
 248
 249         for (fTrack = 0; fTrack < ntracks; fTrack++) {
 250             gAlice->ResetHits();
 251             treeH->GetEntry(fTrack,0);
 252 //
 253 //   Loop over hits
 254             for(AliMUONHit* mHit = (AliMUONHit*)pMUON->FirstHit(-1);
 255                 mHit;
 256                 mHit = (AliMUONHit*)pMUON->NextHit())
 257             {
 258                 fNch = mHit->Chamber()-1;  // chamber number
 259                 if (fNch > AliMUONConstants::NCh()-1) continue;
 260                 iChamber = &(pMUON->Chamber(fNch));
 261
 262 //
 263 // Loop over pad hits
 264                 for (AliMUONPadHit* mPad =
 265                          (AliMUONPadHit*)pMUON->FirstPad(mHit,pMUON->PadHits());
 266                      mPad;
 267                      mPad = (AliMUONPadHit*)pMUON->NextPad(pMUON->PadHits()))
 268                 {
 269                     Int_t cathode  = mPad->Cathode();      // cathode number
 270                     if (cathode != (icat+1)) continue;
 271                     Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
 272 //                  segmentation = iChamber->SegmentationModel(cathode);
 273                     fDigits[0] = mPad->PadX();
 274                     fDigits[1] = mPad->PadY();
 275                     if (!(fHitMap[fNch]->ValidateHit(fDigits[0], fDigits[1]))) continue;
 276                     fDigits[2] = icat;
 277                     fDigits[3] = iqpad;
 278                     fDigits[4] = iqpad;
 279                     if (mHit->Particle() == kMuonPlus ||
 280                         mHit->Particle() == kMuonMinus) {
 281                         fDigits[5] = mPad->HitNumber();
 282                     } else fDigits[5] = -1;
 283
 284                     // build the list of fired pads and update the info
 285
 286                     if (!Exists(mPad)) {
 287                         CreateNew(mPad);
 288                     } else {
 289                         Update(mPad);
 290                     } //  end if pdigit
 291                 } //end loop over clusters
 292             } // hit loop
 293         } // track loop
 294
 295         // open the file with background
 296
 297         if (fMerge) {
 298             fSignal = kFALSE;
 299             ntracks = (Int_t)fTrH1->GetEntries();
 300 //
 301 //   Loop over tracks
 302 //
 303             for (fTrack = 0; fTrack < ntracks; fTrack++) {
 304
 305                 if (fHitsBgr)       fHitsBgr->Clear();
 306                 if (fPadHitsBgr)    fPadHitsBgr->Clear();
 307
 308                 fTrH1->GetEvent(fTrack);
 309 //
 310 //   Loop over hits
 311                 AliMUONHit* mHit;
 312                 for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
 313                 {
 314                     mHit   = (AliMUONHit*) (*fHitsBgr)[i];
 315                     fNch   = mHit->Chamber()-1;  // chamber number
 316                     iChamber = &(pMUON->Chamber(fNch));
 317 //
 318 // Loop over pad hits
 319                     for (AliMUONPadHit* mPad =
 320                              (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHitsBgr);
 321                          mPad;
 322                          mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHitsBgr))
 323                     {
 324                         Int_t cathode  = mPad->Cathode();     // cathode number
 325                         Int_t ipx      = mPad->PadX();        // pad number on X
 326                         Int_t ipy      = mPad->PadY();        // pad number on Y
 327                         Int_t iqpad    = Int_t(mPad->QPad()); // charge per pad
 328                         if (!(fHitMap[fNch]->ValidateHit(ipx, ipy))) continue;
 329
 330                         if (cathode != (icat+1)) continue;
 331                         fDigits[0] = ipx;
 332                         fDigits[1] = ipy;
 333                         fDigits[2] = icat;
 334                         fDigits[3] = iqpad;
 335                         fDigits[4] = 0;
 336                         fDigits[5] = -1;
 337
 338                         // build the list of fired pads and update the info
 339                         if (!Exists(mPad)) {
 340                             CreateNew(mPad);
 341                         } else {
 342                             Update(mPad);
 343                         } //  end if !Exists
 344                     } //end loop over clusters
 345                 } // hit loop
 346             } // track loop
 347
 348             TTree *treeK = gAlice->TreeK();
 349             TFile *file = NULL;
 350
 351             if (treeK) file = treeK->GetCurrentFile();
 352             file->cd();
 353         } // if fMerge
 354
 355         Int_t tracks[kMAXTRACKS];
 356         Int_t charges[kMAXTRACKS];
 357         Int_t nentries = fList->GetEntriesFast();
 358
 359         for (Int_t nent = 0; nent < nentries; nent++) {
 360             AliMUONTransientDigit *address = (AliMUONTransientDigit*)fList->At(nent);
 361             if (address == 0) continue;
 362             Int_t ich = address->Chamber();
 363             Int_t   q = address->Signal();
 364             iChamber = &(pMUON->Chamber(ich));
 365 //
 366 //  Digit Response (noise, threshold, saturation, ...)
 367             AliMUONResponse * response = iChamber->ResponseModel();
 368             q = response->DigitResponse(q,address);
 369
 370             if (!q) continue;
 371
 372             fDigits[0] = address->PadX();
 373             fDigits[1] = address->PadY();
 374             fDigits[2] = address->Cathode();
 375             fDigits[3] = q;
 376             fDigits[4] = address->Physics();
 377             fDigits[5] = address->Hit();
 378
 379             Int_t nptracks = address->GetNTracks();
 380
 381             if (nptracks > kMAXTRACKS) {
 382                 if (fDebug>0)
 383                   printf("\n Attention - nptracks > kMAXTRACKS %d \n", nptracks);
 384                 nptracks = kMAXTRACKS;
 385             }
 386             if (nptracks > 2) {
 387                 if (fDebug>0) {
 388                   printf("Attention - nptracks > 2  %d \n",nptracks);
 389                   printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
 390                 }
 391             }
 392             for (Int_t tr = 0; tr < nptracks; tr++) {
 393                 tracks[tr]   = address->GetTrack(tr);
 394                 charges[tr]  = address->GetCharge(tr);
 395             }      //end loop over list of tracks for one pad
 396             // Sort list of tracks according to charge
 397             if (nptracks > 1) {
 398                 SortTracks(tracks,charges,nptracks);
 399             }
 400             if (nptracks < kMAXTRACKS ) {
 401                 for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
 402                     tracks[i]  = 0;
 403                     charges[i] = 0;
 404                 }
 405             }
 406
 407             // fill digits
 408             pMUON->AddDigits(ich,tracks,charges,fDigits);
 409         }
 410         gAlice->TreeD()->Fill();
 411         pMUON->ResetDigits();
 412         fList->Delete();
 413
 414
 415         for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
 416             if (fHitMap[ii]) {
 417                 delete fHitMap[ii];
 418                 fHitMap[ii] = 0;
 419             }
 420         }
 421     } //end loop over cathodes
 422     delete [] fHitMap;
 423     delete fList;
 424
 425 //  no need to delete ... and it makes a crash also
 426 //    if (fHitsBgr)    fHitsBgr->Delete();
 427 //    if (fPadHitsBgr) fPadHitsBgr->Delete();
 428     // gObjectTable->Print();
 429 }
 430
 431
 432
 433 void AliMUONMerger::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 434 {
 435   //
 436   // Sort the list of tracks contributing to a given digit
 437   // Only the 3 most significant tracks are acctually sorted
 438   //
 439
 440   //
 441   //  Loop over signals, only 3 times
 442   //
 443
 444   Int_t qmax;
 445   Int_t jmax;
 446   Int_t idx[3] = {-2,-2,-2};
 447   Int_t jch[3] = {-2,-2,-2};
 448   Int_t jtr[3] = {-2,-2,-2};
 449   Int_t i,j,imax;
 450
 451   if (ntr<3) imax=ntr;
 452   else imax=3;
 453   for(i=0;i<imax;i++){
 454     qmax=0;
 455     jmax=0;
 456
 457     for(j=0;j<ntr;j++){
 458
 459       if((i == 1 && j == idx[i-1])
 460          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 461
 462       if(charges[j] > qmax) {
 463         qmax = charges[j];
 464         jmax=j;
 465       }
 466     }
 467
 468     if(qmax > 0) {
 469       idx[i]=jmax;
 470       jch[i]=charges[jmax];
 471       jtr[i]=tracks[jmax];
 472     }
 473
 474   }
 475
 476   for(i=0;i<3;i++){
 477     if (jtr[i] == -2) {
 478          charges[i]=0;
 479          tracks[i]=0;
 480     } else {
 481          charges[i]=jch[i];
 482          tracks[i]=jtr[i];
 483     }
 484   }
 485 }
 486
 487
 488