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