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.3  2001/03/05 23:57:44  morsch
  19 Writing of digit tree moved to macro.
  20
  21 Revision 1.2  2001/03/05 08:40:25  morsch
  22 Method SortTracks(..) imported from AliMUON.
  23
  24 Revision 1.1  2001/02/02 14:11:53  morsch
  25 AliMUONMerger prototype to be called by the merge manager.
  26
  27 */
  28
  29 #include <TTree.h>
  30 #include <TVector.h>
  31 #include <TObjArray.h>
  32 #include <TFile.h>
  33 #include <TDirectory.h>
  34
  35
  36 // #include "AliMerger.h"
  37 // #include "AliMergable.h"
  38 #include "AliMUONMerger.h"
  39 #include "AliMUONConstants.h"
  40 #include "AliMUONChamber.h"
  41 #include "AliHitMap.h"
  42 #include "AliMUONHitMapA1.h"
  43 #include "AliMUON.h"
  44 #include "AliMUONHit.h"
  45 #include "AliMUONPadHit.h"
  46 #include "AliMUONDigit.h"
  47 #include "AliMUONTransientDigit.h"
  48 #include "AliRun.h"
  49 #include "AliPDG.h"
  50
  51 ClassImp(AliMUONMerger)
  52
  53 //___________________________________________
  54 AliMUONMerger::AliMUONMerger()
  55 {
  56 // Default constructor
  57     fEvNrSig = 0;
  58     fEvNrBgr = 0;
  59     fMerge   = kDigitize;
  60     fFnBgr   = 0;
  61     fHitMap  = 0;
  62     fList    = 0;
  63     fTrH1    = 0;
  64     fHitsBgr = 0;
  65     fPadHitsBgr = 0;
  66     fHitMap     = 0;
  67     fList       = 0;
  68     fTrList     = 0;
  69     fAddress    = 0;
  70     fBgrFile    = 0;
  71 }
  72
  73 //------------------------------------------------------------------------
  74 AliMUONMerger::~AliMUONMerger()
  75 {
  76 // Destructor
  77     if (fTrH1)       delete fTrH1;
  78     if (fHitsBgr)    delete fHitsBgr;
  79     if (fPadHitsBgr) delete fPadHitsBgr;
  80     if (fHitMap)     delete fHitMap;
  81     if (fList)       delete fList;
  82     if (fTrList)     delete fTrList;
  83     if (fAddress)    delete fAddress;
  84     if (fBgrFile)    delete fBgrFile;
  85 }
  86
  87 //------------------------------------------------------------------------
  88 Bool_t AliMUONMerger::Exists(const AliMUONPadHit *mergable)
  89 {
  90     AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
  91     return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
  92 }
  93
  94 //------------------------------------------------------------------------
  95 void AliMUONMerger::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 AliMUONMerger::CreateNew(AliMUONPadHit *mergable)
 134 {
 135     AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
 136     AliMUONTransientDigit* pdigit;
 137
 138     Int_t ipx      = padhit->PadX();       // pad number on X
 139     Int_t ipy      = padhit->PadY();       // pad number on Y
 140     fList->AddAtAndExpand(
 141         new AliMUONTransientDigit(fNch,fDigits),fCounter);
 142     fHitMap[fNch]->SetHit(ipx, ipy, fCounter);
 143     fCounter++;
 144     pdigit = (AliMUONTransientDigit*)fList->At(fList->GetLast());
 145     // list of tracks
 146     TObjArray *fTrList = (TObjArray*)pdigit->TrackList();
 147     TVector &trinfo    = *((TVector*) (*fAddress)[fCountadr-1]);
 148     fTrList->Add(&trinfo);
 149 }
 150
 151
 152 //------------------------------------------------------------------------
 153 void AliMUONMerger::Init()
 154 {
 155 // Initialisation
 156     // open only once the background file !!
 157     if (fMerge && !fBgrFile) fBgrFile = InitBgr();
 158 }
 159
 160
 161
 162 //------------------------------------------------------------------------
 163 TFile* AliMUONMerger::InitBgr()
 164 {
 165 // Initialise background event
 166     TFile *file = new TFile(fFnBgr);
 167 // add error checking later
 168     printf("\n AliMUONMerger has opened %s file with background event \n", fFnBgr);
 169     fHitsBgr     = new TClonesArray("AliMUONHit",1000);
 170     fPadHitsBgr  = new TClonesArray("AliMUONPadHit",1000);
 171     return file;
 172 }
 173
 174 //------------------------------------------------------------------------
 175 void AliMUONMerger::Digitise()
 176 {
 177
 178     // keep galice.root for signal and name differently the file for
 179     // background when add! otherwise the track info for signal will be lost !
 180
 181     AliMUONChamber*   iChamber;
 182     AliSegmentation*  segmentation;
 183
 184     fList = new TObjArray;
 185     if(!fAddress) fAddress = new TClonesArray("TVector",10000);
 186
 187     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
 188     fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
 189     for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {fHitMap[i] = 0;}
 190     if (fMerge ) {
 191         fBgrFile->cd();
 192         fBgrFile->ls();
 193         //
 194         // Get Hits Tree header from file
 195         //if(fHitsBgr) fHitsBgr->Clear();     // Useless because line 327
 196         //if(fPadHitsBgr) fPadHitsBgr->Clear(); // Useless because line 328
 197         if(fTrH1) delete fTrH1;
 198         fTrH1 = 0;
 199
 200         char treeName[20];
 201         sprintf(treeName,"TreeH%d",fEvNrBgr);
 202         fTrH1 = (TTree*)gDirectory->Get(treeName);
 203         if (!fTrH1) {
 204             printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
 205         }
 206         //
 207         // Set branch addresses
 208         TBranch *branch;
 209         char branchname[20];
 210         sprintf(branchname,"%s",pMUON->GetName());
 211         if (fTrH1 && fHitsBgr) {
 212             branch = fTrH1->GetBranch(branchname);
 213             if (branch) branch->SetAddress(&fHitsBgr);
 214         }
 215         if (fTrH1 && fPadHitsBgr) {
 216             branch = fTrH1->GetBranch("MUONCluster");
 217             if (branch) branch->SetAddress(&fPadHitsBgr);
 218         }
 219     }
 220     //
 221     // loop over cathodes
 222     //
 223     AliHitMap* hm;
 224     fCountadr = 0;
 225     for (int icat = 0; icat < 2; icat++) {
 226         fCounter = 0;
 227         Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
 228         for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
 229             iChamber = &(pMUON->Chamber(i));
 230             if (iChamber->Nsec() == 1 && icat == 1) {
 231                 continue;
 232             } else {
 233                 segmentation = iChamber->SegmentationModel(icat+1);
 234             }
 235             fHitMap[i] = new AliMUONHitMapA1(segmentation, fList);
 236             nmuon[i] = 0;
 237         }
 238
 239 //
 240 //   Loop over tracks
 241 //
 242
 243         TTree *treeH  = gAlice->TreeH();
 244         Int_t ntracks = (Int_t) treeH->GetEntries();
 245         Int_t jj;
 246
 247         Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
 248         for (jj = 0; jj < AliMUONConstants::NCh(); jj++) xhit[jj] = new Float_t[2];
 249         Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
 250         for (jj = 0; jj < AliMUONConstants::NCh(); jj++) yhit[jj] = new Float_t[2];
 251
 252         for (fTrack = 0; fTrack < ntracks; fTrack++) {
 253             gAlice->ResetHits();
 254             treeH->GetEvent(fTrack);
 255 //
 256 //   Loop over hits
 257             for(AliMUONHit* mHit = (AliMUONHit*)pMUON->FirstHit(-1);
 258                 mHit;
 259                 mHit = (AliMUONHit*)pMUON->NextHit())
 260             {
 261                 fNch = mHit->Chamber()-1;  // chamber number
 262                 if (fNch > AliMUONConstants::NCh()-1) continue;
 263                 iChamber = &(pMUON->Chamber(fNch));
 264                 /*
 265                 if (fMerge) {
 266                     if (mHit->Particle() == kMuonPlus ||
 267                         mHit->Particle() == kMuonMinus) {
 268                         xhit[fNch][nmuon[fNch]] = mHit->X();
 269                         yhit[fNch][nmuon[fNch]] = mHit->Y();
 270                         nmuon[fNch]++;
 271                         if (nmuon[fNch] > 2) printf("MUON: nmuon %d\n",nmuon[fNch]);
 272                     }
 273                 }
 274                 */
 275
 276 //
 277 // Loop over pad hits
 278                 for (AliMUONPadHit* mPad =
 279                          (AliMUONPadHit*)pMUON->FirstPad(mHit,pMUON->PadHits());
 280                      mPad;
 281                      mPad = (AliMUONPadHit*)pMUON->NextPad(pMUON->PadHits()))
 282                 {
 283                     Int_t cathode  = mPad->Cathode();      // cathode number
 284                     Int_t ipx      = mPad->PadX();         // pad number on X
 285                     Int_t ipy      = mPad->PadY();         // pad number on Y
 286                     Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
 287                     if (cathode != (icat+1)) continue;
 288
 289                     segmentation = iChamber->SegmentationModel(cathode);
 290
 291                     new((*fAddress)[fCountadr++]) TVector(2);
 292
 293                     TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
 294                     trinfo(0) = (Float_t)fTrack;
 295                     trinfo(1) = (Float_t)iqpad;
 296
 297                     fDigits[0] = ipx;
 298                     fDigits[1] = ipy;
 299                     fDigits[2] = icat;
 300                     fDigits[3] = iqpad;
 301                     fDigits[4] = iqpad;
 302                     if (mHit->Particle() == kMuonPlus ||
 303                         mHit->Particle() == kMuonMinus) {
 304                         fDigits[5] = mPad->HitNumber();
 305                     } else fDigits[5] = -1;
 306
 307                     // build the list of fired pads and update the info
 308
 309                     if (!Exists(mPad)) {
 310                         CreateNew(mPad);
 311                     } else {
 312                         Update(mPad);
 313                     } //  end if pdigit
 314                 } //end loop over clusters
 315             } // hit loop
 316         } // track loop
 317
 318         // open the file with background
 319
 320         if (fMerge) {
 321             ntracks = (Int_t)fTrH1->GetEntries();
 322 //
 323 //   Loop over tracks
 324 //
 325             for (fTrack = 0; fTrack < ntracks; fTrack++) {
 326
 327                 if (fHitsBgr)       fHitsBgr->Clear();
 328                 if (fPadHitsBgr)    fPadHitsBgr->Clear();
 329
 330                 fTrH1->GetEvent(fTrack);
 331 //
 332 //   Loop over hits
 333                 AliMUONHit* mHit;
 334                 for(int i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
 335                 {
 336                     mHit   = (AliMUONHit*) (*fHitsBgr)[i];
 337                     fNch   = mHit->Chamber()-1;  // chamber number
 338                     iChamber = &(pMUON->Chamber(fNch));
 339                     Float_t xbgr = mHit->X();
 340                     Float_t ybgr = mHit->Y();
 341                     Bool_t cond  = kFALSE;
 342                     /*
 343                     for (Int_t imuon = 0; imuon < nmuon[fNch]; imuon++) {
 344                         Float_t dist = (xbgr-xhit[fNch][imuon])*(xbgr-xhit[fNch][imuon])
 345                             +(ybgr-yhit[fNch][imuon])*(ybgr-yhit[fNch][imuon]);
 346                         if (dist < 100.) cond = kTRUE;
 347                     }
 348                     */
 349                     cond  = kTRUE;
 350 //
 351 // Loop over pad hits
 352                     for (AliMUONPadHit* mPad =
 353                              (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHitsBgr);
 354                          mPad;
 355                          mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHitsBgr))
 356                     {
 357                         Int_t cathode  = mPad->Cathode();     // cathode number
 358                         Int_t ipx      = mPad->PadX();        // pad number on X
 359                         Int_t ipy      = mPad->PadY();        // pad number on Y
 360                         Int_t iqpad    = Int_t(mPad->QPad()); // charge per pad
 361
 362                         if (cathode != (icat+1)) continue;
 363                         new((*fAddress)[fCountadr++]) TVector(2);
 364                         TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
 365                         trinfo(0) = kBgTag;  // tag background
 366                         trinfo(1) = kBgTag;
 367
 368                         fDigits[0] = ipx;
 369                         fDigits[1] = ipy;
 370                         fDigits[2] = icat;
 371                         fDigits[3] = iqpad;
 372                         fDigits[4] = 0;
 373                         fDigits[5] = -1;
 374
 375                         // build the list of fired pads and update the info
 376                         if (!Exists(mPad)) {
 377                             CreateNew(mPad);
 378                         } else {
 379                             Update(mPad);
 380                         } //  end if !Exists
 381                     } //end loop over clusters
 382                 } // hit loop
 383             } // track loop
 384
 385             TTree *fAli = gAlice->TreeK();
 386             TFile *file = NULL;
 387
 388             if (fAli) file = fAli->GetCurrentFile();
 389             file->cd();
 390         } // if fMerge
 391         delete [] xhit;
 392         delete [] yhit;
 393
 394         Int_t tracks[10];
 395         Int_t charges[10];
 396         Int_t nentries = fList->GetEntriesFast();
 397
 398         for (Int_t nent = 0; nent < nentries; nent++) {
 399             AliMUONTransientDigit *address = (AliMUONTransientDigit*)fList->At(nent);
 400             if (address == 0) continue;
 401             Int_t ich = address->Chamber();
 402             Int_t   q = address->Signal();
 403             iChamber = &(pMUON->Chamber(ich));
 404 //
 405 //  Digit Response (noise, threshold, saturation, ...)
 406             AliMUONResponse * response = iChamber->ResponseModel();
 407             q = response->DigitResponse(q);
 408
 409             if (!q) continue;
 410
 411             fDigits[0] = address->PadX();
 412             fDigits[1] = address->PadY();
 413             fDigits[2] = address->Cathode();
 414             fDigits[3] = q;
 415             fDigits[4] = address->Physics();
 416             fDigits[5] = address->Hit();
 417
 418             TObjArray* fTrList = (TObjArray*)address->TrackList();
 419             Int_t nptracks = fTrList->GetEntriesFast();
 420
 421             // this was changed to accomodate the real number of tracks
 422
 423             if (nptracks > 10) {
 424                 printf("\n Attention - nptracks > 10 %d \n", nptracks);
 425                 nptracks = 10;
 426             }
 427             if (nptracks > 2) {
 428                 printf("Attention - nptracks > 2  %d \n",nptracks);
 429                 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
 430             }
 431             for (Int_t tr = 0; tr < nptracks; tr++) {
 432                 TVector *ppP = (TVector*)fTrList->At(tr);
 433                 if(!ppP ) printf("ppP - %p\n",ppP);
 434                 TVector &pp  = *ppP;
 435                 tracks[tr]   = Int_t(pp(0));
 436                 charges[tr]  = Int_t(pp(1));
 437             }      //end loop over list of tracks for one pad
 438             // Sort list of tracks according to charge
 439             if (nptracks > 1) {
 440                 SortTracks(tracks,charges,nptracks);
 441             }
 442             if (nptracks < 10 ) {
 443                 for (Int_t i = nptracks; i < 10; i++) {
 444                     tracks[i]  = 0;
 445                     charges[i] = 0;
 446                 }
 447             }
 448
 449             // fill digits
 450             pMUON->AddDigits(ich,tracks,charges,fDigits);
 451             // delete fTrList;
 452         }
 453         gAlice->TreeD()->Fill();
 454         pMUON->ResetDigits();
 455         fList->Delete();
 456
 457
 458         for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
 459             if (fHitMap[ii]) {
 460                 hm=fHitMap[ii];
 461                 delete hm;
 462                 fHitMap[ii] = 0;
 463             }
 464         }
 465         delete [] nmuon;
 466     } //end loop over cathodes
 467     delete [] fHitMap;
 468     delete fList;
 469
 470     if (fAddress)    fAddress->Delete();
 471 //  no need to delete ... and it makes a crash also
 472 //    if (fHitsBgr)    fHitsBgr->Delete();
 473 //    if (fPadHitsBgr) fPadHitsBgr->Delete();
 474     // gObjectTable->Print();
 475 }
 476
 477
 478
 479 void AliMUONMerger::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
 480 {
 481   //
 482   // Sort the list of tracks contributing to a given digit
 483   // Only the 3 most significant tracks are acctually sorted
 484   //
 485
 486   //
 487   //  Loop over signals, only 3 times
 488   //
 489
 490   Int_t qmax;
 491   Int_t jmax;
 492   Int_t idx[3] = {-2,-2,-2};
 493   Int_t jch[3] = {-2,-2,-2};
 494   Int_t jtr[3] = {-2,-2,-2};
 495   Int_t i,j,imax;
 496
 497   if (ntr<3) imax=ntr;
 498   else imax=3;
 499   for(i=0;i<imax;i++){
 500     qmax=0;
 501     jmax=0;
 502
 503     for(j=0;j<ntr;j++){
 504
 505       if((i == 1 && j == idx[i-1])
 506          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
 507
 508       if(charges[j] > qmax) {
 509         qmax = charges[j];
 510         jmax=j;
 511       }
 512     }
 513
 514     if(qmax > 0) {
 515       idx[i]=jmax;
 516       jch[i]=charges[jmax];
 517       jtr[i]=tracks[jmax];
 518     }
 519
 520   }
 521
 522   for(i=0;i<3;i++){
 523     if (jtr[i] == -2) {
 524          charges[i]=0;
 525          tracks[i]=0;
 526     } else {
 527          charges[i]=jch[i];
 528          tracks[i]=jtr[i];
 529     }
 530   }
 531 }
 532
 533
 534