]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ZDC/AliZDCMerger.cxx
Merging and reconstruction code review
[u/mrichter/AliRoot.git] / ZDC / AliZDCMerger.cxx
index e23e3600d0e473590667d1b8767f24bb240bf4cb..a1a6a1454c6f70b8e3ac599fdd8bc936bc05fc21 100644 (file)
 
 ClassImp(AliZDCMerger)
 
+//int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
+
 //____________________________________________________________________________
 AliZDCMerger::AliZDCMerger()
 {
 // Default constructor    
-    fMerge       = kDigitize;
-//   fMerge       = kMerge;
+    //fMerge       = kDigitize -> Only digitization
+    //fMerge       = kMerge    -> Digitization + Merging
+    fMerge       = kMerge;     
     fFnBgr       = 0;
     fBgrFile     = 0;
     fNEvBgr     = 0;
@@ -66,7 +69,7 @@ AliZDCMerger::AliZDCMerger()
     fFnSpecp     = 0;
     fSpecpFile   = 0;
     fNMhits     = 0;
-    fMHits       = 0;
+     
 }
 
 //____________________________________________________________________________
@@ -87,11 +90,10 @@ void AliZDCMerger::InitMerging()
     // Hits tree, impact parameter, num. of spectators n & p 
     //                 in background (full Hijing) event
     Float_t b;
-    Int_t nspecn, nspecp;
-    Background(b, nspecn, nspecp);
+    Background(b, fSpecn, fSpecp);
 
     // Production of nuclear fragments -> num. of FREE spectators n & p
-    Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
+    Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
     
     // Extract from spectators distribution the signal events:
     // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p 
@@ -119,13 +121,18 @@ void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
 //       gAlice = new AliRun("gAlice","Alice test program");
 //      }
 //    }
-
+    
+    gAlice->GetEvent(fNEvBgr);
     AliHeader *header = gAlice->GetHeader();
     AliGenEventHeader* mcHeader = header->GenEventHeader();
     fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
-    fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
-    fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
-    printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
+    Int_t dSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
+    Int_t dSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
+    // Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
+    fSpecn = dSpecn/2;
+    fSpecp = dSpecp/2;
+    printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
+            fNEvBgr,fImpPar,fSpecn,fSpecp);
 }
 
 //____________________________________________________________________________
@@ -143,6 +150,7 @@ TFile* AliZDCMerger::OpenBgrFile()
 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
 {
+    //printf("\n       Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
     Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
     AliZDCFragment *frag = new AliZDCFragment(fImpPar);
     for(j=0; j<=99; j++){
@@ -157,8 +165,11 @@ void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
     Ztot=0;
     Ntot=0;
     frag->AttachNeutrons(zz, nn, Ztot, Ntot);
-    fFreeSpn = fSpecn-Ztot-2*nAlpha;
-    fFreeSpp = fSpecp-Ntot-2*nAlpha;
+    fFreeSpn = fSpecn-Ntot-2*nAlpha;
+    fFreeSpp = fSpecp-Ztot-2*nAlpha;
+    if(fFreeSpn<0) fFreeSpn=0;
+    if(fFreeSpp<0) fFreeSpp=0;
+    //printf("\n                       2*nAlpha = %d, Ztot = %d, Ntot = %d\n",2*nAlpha, Ztot,    Ntot);
     printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
 }
 
@@ -166,7 +177,7 @@ void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
 void AliZDCMerger::Mixing()
 {
     
-//    printf("\n       AliZDCMerger->Mixing\n");
+    //printf("\n       AliZDCMerger->Mixing\n");
        
     // ### Background event Hits ###########################################
     fBgrFile->cd();
@@ -183,7 +194,6 @@ void AliZDCMerger::Mixing()
     fTrSDBgr = 0;    
 
 //    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
-
     // Hits tree
     char treeBgrName[20];
     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
@@ -204,7 +214,7 @@ void AliZDCMerger::Mixing()
     }
     
     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
-//    printf("\n       --- ntracks = %d\n\n", ntracks);
+    //printf("\n       --- ntracks = %d\n\n", ntracks);
     
     Int_t itrack, nhits, ihit, j, sector[2];
     AliZDCHit* zdcHit;
@@ -215,16 +225,9 @@ void AliZDCMerger::Mixing()
 
     // --- Tracks loop
     for(itrack=0; itrack<ntracks; itrack++){
-//       printf("               itrack = %d", itrack);
-//       gAlice->ResetHits();
        fTrHBgr->GetEvent(itrack);
        
-//       for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
-//        zdcHit;
-//        zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
-
        nhits = fHitsBgr->GetEntries();
-//       nhits = ZDCdp->GetNhits();
 //       printf("            nhits = %d \n", nhits);
        for(ihit=0; ihit<nhits; ihit++){    
          zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
@@ -234,7 +237,7 @@ void AliZDCMerger::Mixing()
            MHits[1] = zdcHit->GetXImpact();
            MHits[2] = zdcHit->GetYImpact();
            MHits[3] = zdcHit->GetSFlag();
-           MHits[4] = zdcHit->GetLightPMQ();
+                   MHits[4] = zdcHit->GetLightPMQ();
            MHits[5] = zdcHit->GetLightPMC();
            MHits[6] = zdcHit->GetEnergy();
            MHit = new AliZDCMergedHit(sector, MHits);
@@ -246,11 +249,7 @@ void AliZDCMerger::Mixing()
          }//Hits loop
          
     } // Tracks loop
-//    printf("\n       fNMhits (after bckg) = %d, \n",fNMhits); 
-//    fMHits->Dump();
-//    AliZDCMergedHit *pippo7 = (AliZDCMergedHit*) fMHits->At(7);
-//    pippo7->Dump();
-    
+    //printf(" fNMhits (after bckg) = %d, \n",fNMhits);     
         
     // ### Signal event Hits ###########################################
     // --- Neutrons
@@ -258,7 +257,7 @@ void AliZDCMerger::Mixing()
         
     // --- Protons
     ExtractSignal(2);
-//    printf("\n       fNMhits (after signal) = %d \n",fNMhits); 
+    //printf(" fNMhits (after signal) = %d \n",fNMhits); 
         
 }
 
@@ -273,53 +272,66 @@ void AliZDCMerger::ExtractSignal(Int_t SpecType)
    fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
    fSpecnFile = TFile::Open(fFnSpecn,"R");
    fSpecnFile->cd();
-//   printf("\n        --- ExtractSignal x n: file %s opened\n", fFnSpecn);
+   printf("\n  --- ExtractSignal x n: file %s opened\n", fFnSpecn);
    NumEvents = fFreeSpn;
  }
  else if(SpecType == 2){       // --- Signal for spectator protons
    fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
    fSpecpFile = TFile::Open(fFnSpecp,"R");
    fSpecpFile->cd();
-//   printf("\n        --- ExtractSignal x p: file %s opened\n", fFnSpecp);
+   printf("\n  --- ExtractSignal x p: file %s opened\n", fFnSpecp);
    NumEvents = fFreeSpp;
  }
-// printf("\n          # of free spectator = %d\n", NumEvents);
-// printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
+ //printf("\n          # of free spectators = %d\n", NumEvents);
+ //printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
 
   TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
   Int_t nentries = (Int_t) ZDCSignal->GetEntries();
-//  printf("\n   # entries = %d\n", nentries);
-  Int_t i, j, k;
+  //printf("\n   # entries = %d\n", nentries);
   
   AliZDCMergedHit *MHit; 
   Float_t *entry, HitsSpec[7];
-  Int_t iev, rnd, Volume[2];
-  for(iev=0; iev<NumEvents; iev++){
-     rnd = (Int_t) (1000*gRandom->Rndm());
-//     printf("\n      rnd = %d\n", rnd);
+  Int_t pl, i, j, k, iev=0, rnd[125], Volume[2];
+  for(pl=0;pl<125;pl++){
+     rnd[pl] = 0;
+  }
+  for(pl=0;pl<NumEvents;pl++){
+     rnd[pl] = (Int_t) (9999*gRandom->Rndm());
+     if(rnd[pl] >= 9998) rnd[pl] = 9997;
+     //printf("        rnd[%d] = %d\n",pl,rnd[pl]);     
+  }
+  // Sorting vector in ascending order with C function QSORT 
+  qsort((void*)rnd,NumEvents,sizeof(Int_t),comp);
+  //for(pl=0;pl<NumEvents;pl++){
+  ////printf(" rnd[%d] = %d\n",pl,rnd[pl]);     
+  //}
+  do{
      for(i=0; i<nentries; i++){  
        ZDCSignal->GetEvent(i);
        entry = ZDCSignal->GetArgs();
-       if(entry[0] == rnd){
-//        printf("\n   entry[0] = %f\n", entry[0]);
+       if(entry[0] == rnd[iev]){
           for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
           for(j=0; j<7; j++){
              HitsSpec[j] = entry[j+3];
           }
+          //printf("\n   i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
          MHit = new AliZDCMergedHit(Volume, HitsSpec);
          new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
          delete MHit;
        }
-       else if(entry[0] > rnd) break;
+       else if(entry[0] > rnd[iev]){
+         iev++;
+         continue;
+       }
      }
-  }
+  }while(iev<NumEvents);
   
   if(SpecType ==1){
-//    printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
+    //printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
     fSpecnFile->Close();
   }
   else if(SpecType == 2){
-//    printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
+    //printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
     fSpecpFile->Close();
   }
       
@@ -329,7 +341,7 @@ void AliZDCMerger::ExtractSignal(Int_t SpecType)
 void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
 {
 
-//  printf("\n AliZDCMerger->Digitize()");
+  printf("\n   AliZDCMerger->Digitize()");
 
   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
 //  if(ZDC) printf("\n         Ho trovato lo ZDC!\n");
@@ -344,21 +356,22 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
   
   AliZDCMergedHit *MHit;
   Int_t imhit;
-//  printf("\n    fNMHits = %d\n", fNMhits);
+  printf("        fNMHits = %d\n", fNMhits);
   // Loop over SDigits
   for(imhit=0; imhit<fNMhits; imhit++){
      
      MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
      sector[0] = MHit->GetSector(0);
      sector[1] = MHit->GetSector(1);
+   // tmp -> c'erano i quadranti cannati!
   if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
-//     printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
-//           sector[0], sector[1]);
+     printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
+             sector[0], sector[1]);
        sector[1] = 0;
    }
      lightQ    = Int_t(MHit->GetLightPMQ());
      lightC    = Int_t(MHit->GetLightPMC());
-//     printf("\n    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d", 
+//     printf("    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n", 
 //      imhit,sector[0], sector[1],lightQ, lightC);
      
      if(sector[0] == 1){          //ZN 
@@ -381,7 +394,7 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
   sector[0] = 1; // Detector = ZN
   sector[1] = 0; // Common PM ADC
   digit = Phe2ADCch(1, 0, PMCZN);
-//  printf("\n\n       ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
+  printf("\n\n ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
   PedValue = AddPedestal();
   digit += PedValue;
 //  printf("   PedValue = %d",PedValue);
@@ -390,46 +403,46 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
   for(j=0; j<4; j++){
     sector[1] = j+1; // Towers PM ADCs
     digit = Phe2ADCch(1, j+1, PMQZN[j]);
-//    printf("\n               PMQZN[%d] = %d  phe     ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
+    printf("\n         PMQZN[%d] = %d  phe     ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
     PedValue = AddPedestal();
     digit += PedValue;
 //    printf(" PedValue = %d",PedValue);
     ZDC->AddDigit(sector, digit);
   }
-//  printf("\n");
+  printf("\n");
   
   // Create digits for ZP
   sector[0] = 2; // Detector = ZP
   sector[1] = 0; // Common PM ADC
   digit = Phe2ADCch(2, 0, PMCZP);
-//  printf("\n ZP --- PMCZP = %d       phe     ADCZP = %d ADCch",PMCZP,digit);
+  printf("\n   ZP --- PMCZP = %d       phe     ADCZP = %d ADCch",PMCZP,digit);
   PedValue = AddPedestal();
   digit += PedValue;
-  printf("     PedValue = %d",PedValue);
+//  printf("   PedValue = %d",PedValue);
   ZDC->AddDigit(sector, digit);
   for(j=0; j<4; j++){
     sector[1] = j+1; // Towers PM ADCs
     digit = Phe2ADCch(2, j+1, PMQZP[j]);
-//    printf("\n              PMQZP[%d] = %d   phe     ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
+    printf("\n        PMQZP[%d] = %d   phe     ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
     PedValue = AddPedestal();
     digit += PedValue;
 //    printf(" PedValue = %d",PedValue);
     ZDC->AddDigit(sector, digit);
   }
-//  printf("\n");
+  printf("\n");
   
   // Create digits for ZEM
   sector[0] = 3; 
   sector[1] = 1; // Detector = ZEM1
   digit  = Phe2ADCch(3, 1, PMZEM1);
-//  printf("\n  ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
+  printf("\n   ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
   PedValue = AddPedestal();
   digit += PedValue;
 //  printf("  PedValue = %d\n",PedValue);
   ZDC->AddDigit(sector, digit); 
   sector[1] = 2; // Detector = ZEM2
   digit  = Phe2ADCch(3, 2, PMZEM2);
-//  printf("\n  ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch",PMZEM2,digit);
+  printf("\n   ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch\n",PMZEM2,digit);
   PedValue = AddPedestal();
   digit += PedValue;
 //  printf("  PedValue = %d\n",PedValue);
@@ -446,9 +459,6 @@ Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
   }
   
   Int_t ADCch = 0;
-  // Parameters for conversion of light yield in ADC channels
-  Float_t fPMGain[3][5];      // PM gain
-  Float_t fADCRes;            // ADC conversion factor
   
   Int_t j,i;
   for(i=0; i<3; i++){
@@ -471,7 +481,7 @@ Int_t AliZDCMerger::AddPedestal()
   
   Int_t PedValue;
   Float_t PedMean  = 50.;
-  Float_t PedWidth = 10.;
+  Float_t PedWidth = 5.;
   
   PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);