ZDC/AliZDCMerger.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //                      ZDC event merging class                              //
  20 //                                                                           //
  21 ///////////////////////////////////////////////////////////////////////////////
  22
  23 // --- ROOT system
  24 #include <iostream.h>
  25 #include <TTree.h>
  26 #include <TFile.h>
  27 #include <TDirectory.h>
  28 #include <TNtuple.h>
  29 #include <TSystem.h>
  30
  31 // --- AliRoot header files
  32 #include "AliZDCMerger.h"
  33 #include "AliZDC.h"
  34 #include "AliZDCHit.h"
  35 #include "AliZDCMergedHit.h"
  36 #include "AliZDCDigit.h"
  37 #include "AliZDCFragment.h"
  38 #include "AliRun.h"
  39 #include "AliDetector.h"
  40 #include "AliHeader.h"
  41 #include "AliGenEventHeader.h"
  42 #include "AliGenHijingEventHeader.h"
  43
  44 ClassImp(AliZDCMerger)
  45
  46 //int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
  47
  48 //____________________________________________________________________________
  49 AliZDCMerger::AliZDCMerger()
  50 {
  51 // Default constructor
  52     //fMerge       = kDigitize  -> Only digitization
  53     //fMerge       = kMerge     -> Digitization + Merging
  54     fMerge       = kMerge;
  55     fFnBgr       = 0;
  56     fBgrFile     = 0;
  57     fNEvBgr      = 0;
  58     fTrHBgr      = 0;
  59     fHitsBgr     = 0;
  60     fTrSDBgr     = 0;
  61     fTrDBgr      = 0;
  62     fImpPar      = 0;
  63     fSpecn       = 0;
  64     fSpecp       = 0;
  65     fFreeSpn     = 0;
  66     fFreeSpp     = 0;
  67     fFnSpecn     = 0;
  68     fSpecnFile   = 0;
  69     fFnSpecp     = 0;
  70     fSpecpFile   = 0;
  71     fNMhits      = 0;
  72
  73 }
  74
  75 //____________________________________________________________________________
  76 AliZDCMerger::~AliZDCMerger()
  77 {
  78 // Destructor
  79     if (fBgrFile)    delete fBgrFile;
  80     if (fTrHBgr)     delete fTrHBgr;
  81     if (fTrSDBgr)    delete fTrSDBgr;
  82     if (fHitsBgr)    delete fHitsBgr;
  83     if (fSpecnFile)  delete fSpecnFile;
  84     if (fSpecpFile)  delete fSpecpFile;
  85 }
  86
  87 //____________________________________________________________________________
  88 void AliZDCMerger::InitMerging()
  89 {
  90     // Hits tree, impact parameter, num. of spectators n & p
  91     //          in background (full Hijing) event
  92     Float_t b;
  93     Background(b, fSpecn, fSpecp);
  94
  95     // Production of nuclear fragments -> num. of FREE spectators n & p
  96     Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
  97
  98     // Extract from spectators distribution the signal events:
  99     // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
 100     Mixing();
 101 }
 102
 103 //____________________________________________________________________________
 104 void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
 105 {
 106
 107     // --- Open the background file
 108     if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
 109
 110     // --- Read from the TreeE impact parameter (b),
 111     //     # of spectators n and p (fSpecn, fSpecp)
 112     fBgrFile->cd();
 113
 114 //    // Get AliRun object from file or create it if not on file
 115 //    gAlice = (AliRun*)fBgrFile->Get("gAlice");
 116 //    if (!gAlice) {
 117 //      gAlice = (AliRun*)fBgrFile->Get("gAlice");
 118 //      if (gAlice) printf("AliRun object found on file\n");
 119 //      if (!gAlice) {
 120 //        printf("\n create new gAlice object");
 121 //        gAlice = new AliRun("gAlice","Alice test program");
 122 //      }
 123 //    }
 124
 125     gAlice->GetEvent(fNEvBgr);
 126     AliHeader *header = gAlice->GetHeader();
 127     AliGenEventHeader* mcHeader = header->GenEventHeader();
 128     fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
 129     Int_t dSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
 130     Int_t dSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
 131     // Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
 132     fSpecn = dSpecn/2;
 133     fSpecp = dSpecp/2;
 134     printf("\n  HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
 135             fNEvBgr,fImpPar,fSpecn,fSpecp);
 136 }
 137
 138 //____________________________________________________________________________
 139 TFile* AliZDCMerger::OpenBgrFile()
 140 {
 141     // Initialise background event
 142     TFile *file = new TFile(fFnBgr,"UPDATE");
 143     printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
 144     fHitsBgr = new TClonesArray("AliZDCHit",1000);
 145     fMHits   = new TClonesArray("AliZDCMergedHit",1000);
 146     return file;
 147 }
 148
 149 //____________________________________________________________________________
 150 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
 151                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
 152 {
 153     //printf("\n        Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
 154     Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
 155     AliZDCFragment *frag = new AliZDCFragment(fImpPar);
 156     for(j=0; j<=99; j++){
 157        zz[j] =0;
 158        nn[j] =0;
 159     }
 160
 161     // Fragments generation
 162     frag->GenerateIMF(zz, nAlpha);
 163
 164     // Attach neutrons
 165     Ztot=0;
 166     Ntot=0;
 167     frag->AttachNeutrons(zz, nn, Ztot, Ntot);
 168     fFreeSpn = fSpecn-Ntot-2*nAlpha;
 169     fFreeSpp = fSpecp-Ztot-2*nAlpha;
 170     if(fFreeSpn<0) fFreeSpn=0;
 171     if(fFreeSpp<0) fFreeSpp=0;
 172     //printf("\n                        2*nAlpha = %d, Ztot = %d, Ntot = %d\n",2*nAlpha, Ztot,    Ntot);
 173     printf("\n  Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
 174 }
 175
 176 //____________________________________________________________________________
 177 void AliZDCMerger::Mixing()
 178 {
 179
 180     //printf("\n        AliZDCMerger->Mixing\n");
 181
 182     // ### Background event Hits ###########################################
 183     fBgrFile->cd();
 184 //    fBgrFile->ls();
 185
 186    AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
 187 //    if(ZDC) printf("\n        Ho trovato lo ZDC!\n");
 188
 189     // Hits tree
 190     if(fTrHBgr) delete fTrHBgr;
 191     fTrHBgr = 0;
 192     // SDigits tree
 193     if(fTrSDBgr) delete fTrSDBgr;
 194     fTrSDBgr = 0;
 195
 196 //    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
 197     // Hits tree
 198     char treeBgrName[20];
 199     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
 200     fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
 201     if(!fTrHBgr){
 202       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
 203     }
 204 //    fTrHBgr->Print();
 205
 206     // Branch address
 207     TBranch *branch;
 208     char branchname[20];
 209     sprintf(branchname,"%s",ZDC->GetName());
 210     if(fTrHBgr && fHitsBgr){
 211 //      printf("\n      fTrHBgr!=0 && fHitsBgr!=0\n");
 212       branch = fTrHBgr->GetBranch(branchname);
 213       if(branch) branch->SetAddress(&fHitsBgr);
 214     }
 215
 216     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
 217     //printf("\n        --- ntracks = %d\n\n", ntracks);
 218
 219     Int_t itrack, nhits, ihit, j, sector[2];
 220     AliZDCHit* zdcHit;
 221     AliZDCMergedHit *MHit;
 222     Float_t MHits[7];
 223     TClonesArray &sdigits = *fMHits;    // SDigits TCArray
 224     fNMhits = 0;
 225
 226     // --- Tracks loop
 227     for(itrack=0; itrack<ntracks; itrack++){
 228        fTrHBgr->GetEvent(itrack);
 229
 230        nhits = fHitsBgr->GetEntries();
 231 //       printf("             nhits = %d \n", nhits);
 232        for(ihit=0; ihit<nhits; ihit++){
 233           zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
 234
 235             for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
 236             MHits[0] = zdcHit->GetPrimKinEn();
 237             MHits[1] = zdcHit->GetXImpact();
 238             MHits[2] = zdcHit->GetYImpact();
 239             MHits[3] = zdcHit->GetSFlag();
 240             MHits[4] = zdcHit->GetLightPMQ();
 241             MHits[5] = zdcHit->GetLightPMC();
 242             MHits[6] = zdcHit->GetEnergy();
 243             MHit = new AliZDCMergedHit(sector, MHits);
 244 //          MHit->Print("");
 245             new((*fMHits)[fNMhits]) AliZDCMergedHit(*MHit);
 246             new (sdigits[fNMhits])  AliZDCMergedHit(*MHit);
 247             delete MHit;
 248             fNMhits++;
 249           }//Hits loop
 250
 251     } // Tracks loop
 252     //printf("  fNMhits (after bckg) = %d, \n",fNMhits);
 253
 254     // ### Signal event Hits ###########################################
 255     // --- Neutrons
 256     ExtractSignal(1);
 257
 258     // --- Protons
 259     ExtractSignal(2);
 260     //printf("  fNMhits (after signal) = %d \n",fNMhits);
 261
 262 }
 263
 264 //____________________________________________________________________________
 265 void AliZDCMerger::ExtractSignal(Int_t SpecType)
 266 {
 267
 268 // printf("\n   Entering in Extract Signal\n");
 269
 270  Int_t NumEvents = 0;
 271  if(SpecType == 1){             // --- Signal for spectator neutrons
 272    fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
 273    fSpecnFile = TFile::Open(fFnSpecn,"R");
 274    fSpecnFile->cd();
 275    printf("\n   --- ExtractSignal x n: file %s opened\n", fFnSpecn);
 276    NumEvents = fFreeSpn;
 277  }
 278  else if(SpecType == 2){        // --- Signal for spectator protons
 279    fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
 280    fSpecpFile = TFile::Open(fFnSpecp,"R");
 281    fSpecpFile->cd();
 282    printf("\n   --- ExtractSignal x p: file %s opened\n", fFnSpecp);
 283    NumEvents = fFreeSpp;
 284  }
 285  //printf("\n           # of free spectators = %d\n", NumEvents);
 286  //printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
 287
 288   TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
 289   Int_t nentries = (Int_t) ZDCSignal->GetEntries();
 290   //printf("\n   # entries = %d\n", nentries);
 291
 292   AliZDCMergedHit *MHit;
 293   Float_t *entry, HitsSpec[7];
 294   Int_t pl, i, j, k, iev=0, rnd[125], Volume[2];
 295   for(pl=0;pl<125;pl++){
 296      rnd[pl] = 0;
 297   }
 298   for(pl=0;pl<NumEvents;pl++){
 299      rnd[pl] = (Int_t) (9999*gRandom->Rndm());
 300      if(rnd[pl] >= 9998) rnd[pl] = 9997;
 301      //printf(" rnd[%d] = %d\n",pl,rnd[pl]);
 302   }
 303   // Sorting vector in ascending order with C function QSORT
 304   qsort((void*)rnd,NumEvents,sizeof(Int_t),comp);
 305   //for(pl=0;pl<NumEvents;pl++){
 306   ////printf("  rnd[%d] = %d\n",pl,rnd[pl]);
 307   //}
 308   do{
 309      for(i=0; i<nentries; i++){
 310         ZDCSignal->GetEvent(i);
 311         entry = ZDCSignal->GetArgs();
 312         if(entry[0] == rnd[iev]){
 313           for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
 314           for(j=0; j<7; j++){
 315              HitsSpec[j] = entry[j+3];
 316           }
 317           //printf("\n   i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
 318           MHit = new AliZDCMergedHit(Volume, HitsSpec);
 319           new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
 320           delete MHit;
 321         }
 322         else if(entry[0] > rnd[iev]){
 323           iev++;
 324           continue;
 325         }
 326      }
 327   }while(iev<NumEvents);
 328
 329   if(SpecType ==1){
 330     //printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
 331     fSpecnFile->Close();
 332   }
 333   else if(SpecType == 2){
 334     //printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
 335     fSpecpFile->Close();
 336   }
 337
 338 }
 339
 340 //____________________________________________________________________________
 341 void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
 342 {
 343
 344   printf("\n    AliZDCMerger->Digitize()");
 345
 346   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
 347 //  if(ZDC) printf("\n  Ho trovato lo ZDC!\n");
 348
 349   Int_t lightQ, lightC, sector[2], digit;
 350   Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
 351   Int_t i;
 352   for(i=0; i<4; i++){
 353      PMQZN[i] = 0;
 354      PMQZP[i] = 0;
 355   }
 356
 357   AliZDCMergedHit *MHit;
 358   Int_t imhit;
 359   printf("         fNMHits = %d\n", fNMhits);
 360   // Loop over SDigits
 361   for(imhit=0; imhit<fNMhits; imhit++){
 362
 363      MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
 364      sector[0] = MHit->GetSector(0);
 365      sector[1] = MHit->GetSector(1);
 366    // tmp -> c'erano i quadranti cannati!
 367   if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
 368      printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
 369               sector[0], sector[1]);
 370        sector[1] = 0;
 371    }
 372      lightQ    = Int_t(MHit->GetLightPMQ());
 373      lightC    = Int_t(MHit->GetLightPMC());
 374 //     printf("    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n",
 375 //      imhit,sector[0], sector[1],lightQ, lightC);
 376
 377      if(sector[0] == 1){           //ZN
 378        PMCZN = PMCZN + lightC;
 379        PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
 380      }
 381      else if(sector[0] == 2){      //ZP
 382        PMCZP = PMCZP + lightC;
 383        PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
 384      }
 385      else if(sector[0] == 3){      //ZEM
 386        if(sector[1] ==1) PMZEM1 = PMZEM1 + lightC;
 387        else              PMZEM2 = PMZEM2 + lightQ;
 388      }
 389   } // SDigits loop
 390
 391   // ### Digits creation ###############################################
 392   // Create digits for ZN
 393   Int_t PedValue;
 394   sector[0] = 1; // Detector = ZN
 395   sector[1] = 0; // Common PM ADC
 396   digit = Phe2ADCch(1, 0, PMCZN);
 397   printf("\n\n  ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
 398   PedValue = AddPedestal();
 399   digit += PedValue;
 400 //  printf("    PedValue = %d",PedValue);
 401   ZDC->AddDigit(sector, digit);
 402   Int_t j;
 403   for(j=0; j<4; j++){
 404     sector[1] = j+1; // Towers PM ADCs
 405     digit = Phe2ADCch(1, j+1, PMQZN[j]);
 406     printf("\n          PMQZN[%d] = %d  phe     ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
 407     PedValue = AddPedestal();
 408     digit += PedValue;
 409 //    printf("  PedValue = %d",PedValue);
 410     ZDC->AddDigit(sector, digit);
 411   }
 412   printf("\n");
 413
 414   // Create digits for ZP
 415   sector[0] = 2; // Detector = ZP
 416   sector[1] = 0; // Common PM ADC
 417   digit = Phe2ADCch(2, 0, PMCZP);
 418   printf("\n    ZP --- PMCZP = %d       phe     ADCZP = %d ADCch",PMCZP,digit);
 419   PedValue = AddPedestal();
 420   digit += PedValue;
 421 //  printf("    PedValue = %d",PedValue);
 422   ZDC->AddDigit(sector, digit);
 423   for(j=0; j<4; j++){
 424     sector[1] = j+1; // Towers PM ADCs
 425     digit = Phe2ADCch(2, j+1, PMQZP[j]);
 426     printf("\n         PMQZP[%d] = %d   phe     ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
 427     PedValue = AddPedestal();
 428     digit += PedValue;
 429 //    printf("  PedValue = %d",PedValue);
 430     ZDC->AddDigit(sector, digit);
 431   }
 432   printf("\n");
 433
 434   // Create digits for ZEM
 435   sector[0] = 3;
 436   sector[1] = 1; // Detector = ZEM1
 437   digit  = Phe2ADCch(3, 1, PMZEM1);
 438   printf("\n    ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
 439   PedValue = AddPedestal();
 440   digit += PedValue;
 441 //  printf("  PedValue = %d\n",PedValue);
 442   ZDC->AddDigit(sector, digit);
 443   sector[1] = 2; // Detector = ZEM2
 444   digit  = Phe2ADCch(3, 2, PMZEM2);
 445   printf("\n    ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch\n",PMZEM2,digit);
 446   PedValue = AddPedestal();
 447   digit += PedValue;
 448 //  printf("  PedValue = %d\n",PedValue);
 449   ZDC->AddDigit(sector, digit);
 450 }
 451
 452 //_____________________________________________________________________________
 453 Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
 454 {
 455   // Evaluation of the ADC channel corresponding to the light yield Light
 456
 457   if(gAlice->GetDebug() > 0){
 458 //    printf("\n  Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
 459   }
 460
 461   Int_t ADCch = 0;
 462
 463   Int_t j,i;
 464   for(i=0; i<3; i++){
 465      for(j=0; j<5; j++){
 466         fPMGain[i][j]   = 100000.;
 467      }
 468   }
 469   fADCRes   = 0.00000064; // ADC Resolution: 250 fC/ADCch
 470
 471   ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
 472
 473   return ADCch;
 474 }
 475
 476 //_____________________________________________________________________________
 477 Int_t AliZDCMerger::AddPedestal()
 478 {
 479   // --- Pedestal value -> extracted from a gaussian distribution
 480   // obtained from the beam test on the ZEM prototype (Aug. 2000)
 481
 482   Int_t PedValue;
 483   Float_t PedMean  = 50.;
 484   Float_t PedWidth = 5.;
 485
 486   PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);
 487
 488   return PedValue;
 489 }