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