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