Merge branch 'master' of https://git.cern.ch/reps/AliRoot

diff --git a/EMCAL/AliEMCALDigitizer.cxx b/EMCAL/AliEMCALDigitizer.cxx
index f7f72858d730bdef5f4364327f38372db32cdbb6..2524e99a70d873cc0f0fe6ca9dcdbdd749efd67f 100644 (file)
--- a/EMCAL/AliEMCALDigitizer.cxx
+++ b/EMCAL/AliEMCALDigitizer.cxx
@@ -268,364 +268,373 @@ void AliEMCALDigitizer::Digitize(Int_t event)
   // must call the Calibrate() method of the SDigitizer to convert it
   // back to an energy in GeV before adding it to the Digit
   //
-  static int nEMC=0; //max number of digits possible
+  
   AliRunLoader *rl = AliRunLoader::Instance();
   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
   
-  if(emcalLoader)
+  if(!emcalLoader)
   {
-    Int_t readEvent = event ;
-    // fDigInput is data member from AliDigitizer
-    if (fDigInput)
+    AliFatal("EMCALLoader is NULL!") ;
+    return; // not needed but in case coverity complains ...
+  }
+  
+  Int_t readEvent = event ;
+  if (fDigInput) // fDigInput is data member from AliDigitizer
     readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
-    AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
-                    readEvent, GetTitle(), fEventFolderName.Data())) ;
-    rl->GetEvent(readEvent);
+  AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
+                  readEvent, GetTitle(), fEventFolderName.Data())) ;
+  rl->GetEvent(readEvent);
+  
+  TClonesArray * digits = emcalLoader->Digits() ;
+  digits->Delete() ;  //correct way to clear array when memory is
+  //allocated by objects stored in array
+  
+  // Load Geometry
+  if (!rl->GetAliRun())
+  {
+    AliFatal("Could not get AliRun from runLoader");
+    return; // not needed but in case coverity complains ...
+  }
+  
+  AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
+  AliEMCALGeometry *geom = emcal->GetGeometry();
+  static int nEMC = geom->GetNCells();//max number of digits possible
+  AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
+  
+  digits->Expand(nEMC) ;
+  
+  // RS create a digitizer on the fly
+  if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
+  fSDigitizer->SetEventRange(0, -1) ;
+  
+  //-------------------------------------------------------
+  //take all the inputs to add together and load the SDigits
+  TObjArray * sdigArray = new TObjArray(fInput) ;
+  sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
+  
+  Int_t i ;
+  Int_t embed = kFALSE;
+  for(i = 1 ; i < fInput ; i++)
+  {
+    TString tempo(fEventNames[i]) ;
+    tempo += i ;
+    
+    AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
+    
+    if (!rl2)
+      rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
     
-    TClonesArray * digits = emcalLoader->Digits() ;
-    digits->Delete() ;  //correct way to clear array when memory is
-    //allocated by objects stored in array
+    if(!rl2)
+    {
+      AliFatal("Run Loader of second event not available!");
+      return; // not needed but in case coverity complains ...
+    }
     
-    // Load Geometry
-    AliEMCALGeometry *geom = 0;
-    if (!rl->GetAliRun())
+    if (fDigInput)
+      readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
+    
+    Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
+    
+    rl2->LoadSDigits();
+    //     rl2->LoadDigits();
+    rl2->GetEvent(readEvent);
+    
+    if(!rl2->GetDetectorLoader("EMCAL"))
     {
-      AliFatal("Could not get AliRun from runLoader");
+      AliFatal("Loader of second input not found");
+      return; // not needed but in case coverity complains ...
     }
-    else{
-      AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
-      geom = emcal->GetGeometry();
-      nEMC = geom->GetNCells();
-      AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
+    
+    AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
+    
+    if(!emcalLoader2)
+    {
+      AliFatal("EMCAL Loader of second event not available!");
+      return; // not needed but in case coverity complains ...
     }
+
+    //Merge 2 simulated sdigits
+    if(!emcalLoader2->SDigits()) continue;
+    
+    TClonesArray* sdigits2 = emcalLoader2->SDigits();
+    sdigArray->AddAt(sdigits2, i) ;
     
-    Int_t absID = -1 ;
+    // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
+    // do not smear energy of embedded, do not add noise to any sdigits
+    if( sdigits2->GetEntriesFast() <= 0 ) continue;
     
-    digits->Expand(nEMC) ;
+    //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
+    AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
+    if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
     
-    // RS create a digitizer on the fly
-    if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
-    fSDigitizer->SetEventRange(0, -1) ;
-    //
-    //take all the inputs to add together and load the SDigits
-    TObjArray * sdigArray = new TObjArray(fInput) ;
-    sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
-    Int_t i ;
-    Int_t embed = kFALSE;
-    for(i = 1 ; i < fInput ; i++)
+  }// input loop
+  
+  //--------------------------------
+  //Find the first tower with signal
+  Int_t nextSig = nEMC + 1 ;
+  TClonesArray * sdigits ;
+  for(i = 0 ; i < fInput ; i++)
+  {
+    if(i > 0 && embed) continue;
+    
+    sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
+    if(!sdigits)
     {
-      TString tempo(fEventNames[i]) ;
-      tempo += i ;
-      
-      AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
-      
-      if (rl2==0)
-      rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
-      
-      if (fDigInput)
-      readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
-      Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
-      rl2->LoadSDigits();
-      //     rl2->LoadDigits();
-      rl2->GetEvent(readEvent);
-      if(rl2->GetDetectorLoader("EMCAL"))
-      {
-        AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
-        
-        if(emcalLoader2)
-        {
-          //Merge 2 simulated sdigits
-          if(emcalLoader2->SDigits()){
-            TClonesArray* sdigits2 = emcalLoader2->SDigits();
-            sdigArray->AddAt(sdigits2, i) ;
-            // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
-            // do not smear energy of embedded, do not add noise to any sdigits
-            if(sdigits2->GetEntriesFast()>0)
-            {
-              //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
-              AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
-              if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
-            }
-          }
-          
-        }//loader2
-        else  AliFatal("EMCAL Loader of second event not available!");
-        
-      }
-      else Fatal("Digitize", "Loader of second input not found");
-    }// input loop
+      AliDebug(1,"Null sdigit pointer");
+      continue;
+    }
     
-    //Find the first tower with signal
-    Int_t nextSig = nEMC + 1 ;
-    TClonesArray * sdigits ;
-    for(i = 0 ; i < fInput ; i++){
-      if(i > 0 && embed) continue;
-      sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
-      if(sdigits)
-      {
-        if (sdigits->GetEntriesFast() )
-        {
-          AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
-          if(sd)
-          {
-            Int_t curNext = sd->GetId() ;
-            if(curNext < nextSig)
-            nextSig = curNext ;
-            AliDebug(1,Form("input %i : #sdigits %i \n",
-                            i, sdigits->GetEntriesFast()));
-            
-          }//First entry is not NULL
-          else
-          {
-            AliDebug(1, "NULL sdigit pointer");
-            continue;
-          }
-        }//There is at least one entry
-        else
-        {
-          AliDebug(1, "NULL sdigits array");
-          continue;
-        }
-      }// SDigits array exists
-      else
-      {
-        AliDebug(1,"Null sdigit pointer");
-        continue;
-      }
-    }// input loop
-    AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
+    if (sdigits->GetEntriesFast() <= 0 )
+    {
+      AliDebug(1, "No sdigits entries");
+      continue;
+    }
     
-    TArrayI index(fInput) ;
-    index.Reset() ;  //Set all indexes to zero
+    AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
+    if(!sd)
+    {
+      AliDebug(1, "NULL sdigit pointer");
+      continue;
+    }
     
-    AliEMCALDigit * digit ;
-    AliEMCALDigit * curSDigit ;
+    Int_t curNext = sd->GetId() ;
+    if(curNext < nextSig)
+    nextSig = curNext ;
+    AliDebug(1,Form("input %i : #sdigits %i \n",i, sdigits->GetEntriesFast()));
     
-    //Put Noise contribution, smear time and energy
-    Float_t timeResolution = 0;
-    for(absID = 0; absID < nEMC; absID++)
-    { // Nov 30, 2006 by PAI; was from 1 to nEMC
-      
-      if (IsDead(absID)) continue; // Don't instantiate dead digits
-      
-      Float_t energy = 0 ;
-      
-      // amplitude set to zero, noise will be added later
-      Float_t noiseTime = 0.;
-      if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
-      new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
-      //look if we have to add signal?
-      digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
+  }// input loop
+  
+  AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
+  
+  TArrayI index(fInput) ;
+  index.Reset() ;  //Set all indexes to zero
+  
+  AliEMCALDigit * digit ;
+  AliEMCALDigit * curSDigit ;
+  
+  //---------------------------------------------
+  //Put Noise contribution, smear time and energy
+  Float_t timeResolution = 0;
+  Int_t absID = -1 ;
+  for(absID = 0; absID < nEMC; absID++)
+  { // Nov 30, 2006 by PAI; was from 1 to nEMC
+    
+    if (IsDead(absID)) continue; // Don't instantiate dead digits
+    
+    Float_t energy = 0 ;
+    
+    // amplitude set to zero, noise will be added later
+    Float_t noiseTime = 0.;
+    if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
+    new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
+    //look if we have to add signal?
+    digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
+    
+    if (!digit)
+    {
+      AliDebug(1,"Digit pointer is null");
+      continue;
+    }
+    
+    if(absID==nextSig)
+    {
+      // Calculate time as time of the largest digit
+      Float_t time = digit->GetTime() ;
+      Float_t aTime= digit->GetAmplitude() ;
       
-      if (digit)
+      // loop over input
+      Int_t nInputs = fInput;
+      if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
+      for(i = 0; i< nInputs ; i++)
       {
-        if(absID==nextSig)
+        //loop over (possible) merge sources
+        TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
+        if(sdtclarr)
         {
-          // Calculate time as time of the largest digit
-          Float_t time = digit->GetTime() ;
-          Float_t aTime= digit->GetAmplitude() ;
-          
-          // loop over input
-          Int_t nInputs = fInput;
-          if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
-          for(i = 0; i< nInputs ; i++)
-          {  //loop over (possible) merge sources
-            TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
-            if(sdtclarr) {
-              Int_t sDigitEntries = sdtclarr->GetEntriesFast();
-              
-              if(sDigitEntries > index[i] )
-              curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
-              else
-              curSDigit = 0 ;
-              
-              //May be several digits will contribute from the same input
-              while(curSDigit && (curSDigit->GetId() == absID))
-              {
-                //Shift primary to separate primaries belonging different inputs
-                Int_t primaryoffset ;
-                if(fDigInput)
-                primaryoffset = fDigInput->GetMask(i) ;
-                else
-                primaryoffset = i ;
-                curSDigit->ShiftPrimary(primaryoffset) ;
-                
-                if(curSDigit->GetAmplitude()>aTime)
-                {
-                  aTime = curSDigit->GetAmplitude();
-                  time  = curSDigit->GetTime();
-                }
-                
-                *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
-                
-                index[i]++ ;
-                
-                if( sDigitEntries > index[i] )
-                curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
-                else
-                curSDigit = 0 ;
-              }// while
-            }// source exists
-          }// loop over merging sources
+          Int_t sDigitEntries = sdtclarr->GetEntriesFast();
           
+          if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
+          else                          curSDigit = 0 ;
           
-          //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
-          energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
-          
-          // add fluctuations for photo-electron creation
-          // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
-          Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
-          energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
-          
-          //calculate and set time
-          digit->SetTime(time) ;
-          
-          //Find next signal module
-          nextSig = nEMC + 1 ;
-          for(i = 0 ; i < nInputs ; i++){
-            sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
+          //May be several digits will contribute from the same input
+          while(curSDigit && (curSDigit->GetId() == absID))
+          {
+            //Shift primary to separate primaries belonging different inputs
+            Int_t primaryoffset = i ;
+            if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
+            curSDigit->ShiftPrimary(primaryoffset) ;
             
-            if(sdigits){
-              Int_t curNext = nextSig ;
-              if(sdigits->GetEntriesFast() > index[i])
-              {
-                AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
-                if(tmpdigit)
-                {
-                  curNext = tmpdigit->GetId() ;
-                }
-              }
-              if(curNext < nextSig) nextSig = curNext ;
-            }// sdigits exist
-          } // input loop
-          
-        }//absID==nextSig
-        
-        // add the noise now, no need for embedded events with real data
-        if(!embed)
-        energy += gRandom->Gaus(0., fPinNoise) ;
-        
-        
-        // JLK 26-June-2008
-        //Now digitize the energy according to the fSDigitizer method,
-        //which merely converts the energy to an integer.  Later we will
-        //check that the stored value matches our allowed dynamic ranges
-        digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
-        
-        //Set time resolution with final energy
-        timeResolution = GetTimeResolution(energy);
-        digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
-        AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
-                         absID, energy, nextSig));
-        //Add delay to time
-        digit->SetTime(digit->GetTime()+fTimeDelay) ;
-        
-      }// Digit pointer exists
-      else AliDebug(1,"Digit pointer is null");
+            if(curSDigit->GetAmplitude()>aTime)
+            {
+              aTime = curSDigit->GetAmplitude();
+              time  = curSDigit->GetTime();
+            }
+            
+            *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
+            
+            index[i]++ ;
+            
+            if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
+            else                           curSDigit = 0 ;
+          }// while
+        }// source exists
+      }// loop over merging sources
       
-    } // for(absID = 0; absID < nEMC; absID++)
-    
-    //Embed simulated amplitude (and time?) to real data digits
-    if(embed)
-    {
-      for(Int_t input = 1; input<fInput; input++)
+      //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
+      energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
+      
+      // add fluctuations for photo-electron creation
+      // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
+      Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
+      energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
+      
+      //calculate and set time
+      digit->SetTime(time) ;
+      
+      //Find next signal module
+      nextSig = nEMC + 1 ;
+      for(i = 0 ; i < nInputs ; i++)
       {
-        TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
-        if(!realDigits)
-        {
-          AliDebug(1,"No sdigits to merge\n");
-          continue;
-        }
+        sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
         
-        for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
+        if(sdigits)
         {
-          AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
-          if(digit2)
+          Int_t curNext = nextSig ;
+          if(sdigits->GetEntriesFast() > index[i])
           {
-            digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
-            if(digit)
+            AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
+            if(tmpdigit)
             {
-              // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
-              // and multiply to get a big int.
-              Float_t time2 = digit2->GetTime();
-              Float_t e2    = digit2->GetAmplitude();
-              CalibrateADCTime(e2,time2,digit2->GetId());
-              Float_t amp2 = fSDigitizer->Digitize(e2);
-              
-              Float_t e0    = digit ->GetAmplitude();
-              if(e0>0.01)
-              AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
-                              digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
-                              digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
-              
-              // Sum amplitudes, change digit amplitude
-              digit->SetAmplitude( digit->GetAmplitude() + amp2);
-              // Rough assumption, assign time of the largest of the 2 digits,
-              // data or signal to the final digit.
-              if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
-              //Mark the new digit as embedded
-              digit->SetType(AliEMCALDigit::kEmbedded);
-              
-              if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
-              AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
-                              digit->GetId(), digit->GetAmplitude(), digit->GetType()));
-            }//digit2
-          }//digit2
-        }//loop digit 2
-      }//input loop
-    }//embed
-    
-    //JLK is it better to call Clear() here?
-    delete sdigArray ; //We should not delete its contents
+              curNext = tmpdigit->GetId() ;
+            }
+          }
+          if(curNext < nextSig) nextSig = curNext ;
+        }// sdigits exist
+      } // input loop
+      
+    }//absID==nextSig
     
-    //remove digits below thresholds
-    // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
-    // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
-    // before merge in the same loop real data digits if available
-    Float_t energy = 0;
-    Float_t time   = 0;
-    for(i = 0 ; i < nEMC ; i++)
-    {
-      digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
-      if(digit)
-      {
-        //Then get the energy in GeV units.
-        energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
-        //Then digitize using the calibration constants of the ocdb
-        Float_t ampADC = energy;
-        DigitizeEnergyTime(ampADC, time, digit->GetId())  ;      
-        if(ampADC < fDigitThreshold)
-        digits->RemoveAt(i) ;
-        
-      }// digit exists
-    } // digit loop
+    // add the noise now, no need for embedded events with real data
+    if(!embed)
+      energy += gRandom->Gaus(0., fPinNoise) ;
     
-    digits->Compress() ;  
+    // JLK 26-June-2008
+    //Now digitize the energy according to the fSDigitizer method,
+    //which merely converts the energy to an integer.  Later we will
+    //check that the stored value matches our allowed dynamic ranges
+    digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
     
-    Int_t ndigits = digits->GetEntriesFast() ; 
+    //Set time resolution with final energy
+    timeResolution = GetTimeResolution(energy);
+    digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
+    AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
+                     absID, energy, nextSig));
+    //Add delay to time
+    digit->SetTime(digit->GetTime()+fTimeDelay) ;
     
-    //JLK 26-June-2008
-    //After we have done the summing and digitizing to create the
-    //digits, now we want to calibrate the resulting amplitude to match
-    //the dynamic range of our real data.  
-    for (i = 0 ; i < ndigits ; i++)
+  } // for(absID = 0; absID < nEMC; absID++)
+  
+  //---------------------------------------------------------
+  //Embed simulated amplitude (and time?) to real data digits
+  if(embed)
+  {
+    for(Int_t input = 1; input<fInput; input++)
     {
-      digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
-      if(digit)
+      TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
+      if(!realDigits)
       {
-        digit->SetIndexInList(i) ; 
-        energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
-        time   = digit->GetTime();
-        Float_t ampADC = energy;
-        DigitizeEnergyTime(ampADC, time, digit->GetId());
-        digit->SetAmplitude(ampADC) ;
-        digit->SetTime(time);
-        // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
-      }// digit exists
-    }//Digit loop
-    
-  }
-  else AliFatal("EMCALLoader is NULL!") ;
+        AliDebug(1,"No sdigits to merge\n");
+        continue;
+      }
+      
+      for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
+      {
+        AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
+        if(digit2)
+        {
+          digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
+          if(digit)
+          {
+            // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
+            // and multiply to get a big int.
+            Float_t time2 = digit2->GetTime();
+            Float_t e2    = digit2->GetAmplitude();
+            CalibrateADCTime(e2,time2,digit2->GetId());
+            Float_t amp2 = fSDigitizer->Digitize(e2);
+            
+            Float_t e0    = digit ->GetAmplitude();
+            if(e0>0.01)
+            AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
+                            digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
+                            digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
+            
+            // Sum amplitudes, change digit amplitude
+            digit->SetAmplitude( digit->GetAmplitude() + amp2);
+            // Rough assumption, assign time of the largest of the 2 digits,
+            // data or signal to the final digit.
+            if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
+            //Mark the new digit as embedded
+            digit->SetType(AliEMCALDigit::kEmbedded);
+            
+            if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
+            AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
+                            digit->GetId(), digit->GetAmplitude(), digit->GetType()));
+          }//digit2
+        }//digit2
+      }//loop digit 2
+    }//input loop
+  }//embed
+  
+  //JLK is it better to call Clear() here?
+  delete sdigArray ; //We should not delete its contents
+  
+  //------------------------------
+  //remove digits below thresholds
+  // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
+  // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
+  // before merge in the same loop real data digits if available
+  Float_t energy = 0;
+  Float_t time   = 0;
+  for(i = 0 ; i < nEMC ; i++)
+  {
+    digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
+    if(digit)
+    {
+      //Then get the energy in GeV units.
+      energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
+      //Then digitize using the calibration constants of the ocdb
+      Float_t ampADC = energy;
+      DigitizeEnergyTime(ampADC, time, digit->GetId())  ;
+      if(ampADC < fDigitThreshold)
+      digits->RemoveAt(i) ;
+      
+    }// digit exists
+  } // digit loop
+  
+  digits->Compress() ;
+  
+  Int_t ndigits = digits->GetEntriesFast() ;
+  
+  //---------------------------------------------------------------
+  //JLK 26-June-2008
+  //After we have done the summing and digitizing to create the
+  //digits, now we want to calibrate the resulting amplitude to match
+  //the dynamic range of our real data.
+  for (i = 0 ; i < ndigits ; i++)
+  {
+    digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
+    if(digit)
+    {
+      digit->SetIndexInList(i) ;
+      energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
+      time   = digit->GetTime();
+      Float_t ampADC = energy;
+      DigitizeEnergyTime(ampADC, time, digit->GetId());
+      digit->SetAmplitude(ampADC) ;
+      digit->SetTime(time);
+      // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
+    }// digit exists
+  }//Digit loop
   
 }
 

diff --git a/EMCAL/macros/Config.C b/EMCAL/macros/Config.C
old mode 100755 (executable)
new mode 100644 (file)
index 216e768..80e8fc2
--- a/EMCAL/macros/Config.C
+++ b/EMCAL/macros/Config.C
@@ -154,8 +154,10 @@ void Config()
     gener->SetPhiRange(80.0,120.0);
   else if(year == 2011)
     gener->SetPhiRange(80.0,180.0);
-  else
+  else if(year == 2012 || year == 2013)
     gener->SetPhiRange(80.0,190.0);
+  else
+    gener->SetPhiRange(80.0,330.0); // Include DCal
   
   gener->SetThetaRange(EtaToTheta(0.7), EtaToTheta(-0.7));
   
@@ -378,14 +380,15 @@ void Config()
   if (iEMCAL)
   {
     //=================== EMCAL parameters ============================
+    AliEMCAL *EMCAL = 0;
     if      (year == 2010)  // d phi = 40 degrees
-      AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1",    checkGeoAndRun);
+      EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1",    checkGeoAndRun);
     else if (year == 2011)  // d phi = 100 degrees
-      AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETEV1",     checkGeoAndRun);
-    else if (year > 2011)   // d phi = 110 degrees
-      AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE12SMV1", checkGeoAndRun);
-    else // Old configuration with 110 degrees but not perfect geometry
-      AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE",       checkGeoAndRun);
+      EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETEV1",     checkGeoAndRun);
+    else if (year == 2012 || year == 2013)   // d phi = 107 degrees
+      EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE12SMV1", checkGeoAndRun);
+    else
+      EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE12SMV1_DCAL_8SM", checkGeoAndRun); // EMCAL+DCAL dphi = 107 (EMCAL) + 33 (gap) + 67 (DCAL)
   }
   
   if (iACORDE)

diff --git a/EMCAL/macros/TestEMCALReconstruction.C b/EMCAL/macros/TestEMCALReconstruction.C
index 32e0e45207c3eb9a4f32a67adda27d72cb58b6cb..e9aabf961ea187712465e27f141db00ecaeb1e12 100644 (file)
--- a/EMCAL/macros/TestEMCALReconstruction.C
+++ b/EMCAL/macros/TestEMCALReconstruction.C
@@ -21,11 +21,19 @@ void TestEMCALReconstruction(Int_t nev =-1) {
   rec.SetRunLocalReconstruction("EMCAL");  //only do emcal
   rec.SetFillESD("EMCAL");
   rec.SetEventRange(0,nev);
-  //rec.SetInput("raw.root");
-  //rec.SetRunQA(":");
+  rec.SetRunQA(":");
 
+  // Decomment this line in case of real data,
+  // add the proper name of the file
+  //rec.SetInput("raw.root");
+  
   //OCDB settings
   rec.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  
+  // Decommetn this line in case of anchored MC runs or data,
+  // with the appropriate year
+  //rec.SetDefaultStorage("alien://Folder=/alice/data/2011/OCDB");
+
   rec.SetSpecificStorage("GRP/GRP/Data",
                          Form("local://%s",gSystem->pwd()));
 

diff --git a/EMCAL/macros/TestEMCALSDigit.C b/EMCAL/macros/TestEMCALSDigit.C
new file mode 100644 (file)
index 0000000..64261db
--- /dev/null
+++ b/EMCAL/macros/TestEMCALSDigit.C
@@ -0,0 +1,93 @@
+// Test Macro, shows how to load Digits and Geometry, and how can we get 
+// some of the parameters and variables.
+// Author: Gustavo Conesa
+
+void TestEMCALSDigit()
+{
+  
+  // Getting EMCAL Detector and Geometry.
+  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
+  
+  if (rl == 0x0)
+    cout<<"Can not instatiate the Run Loader"<<endl;
+  
+  rl->LoadgAlice();//Needed to get geometry
+  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
+    (rl->GetDetectorLoader("EMCAL"));
+  
+  TGeoManager::Import("geometry.root");
+  
+  AliRun * alirun   = rl->GetAliRun(); // Needed to get Geometry
+  AliEMCALGeometry * geom ;
+  if(alirun){
+    AliEMCAL * emcal  = (AliEMCAL*)alirun->GetDetector("EMCAL");
+    geom = emcal->GetGeometry();
+  }
+  
+  if (geom == 0) cout<<"Did not get geometry from EMCALLoader"<<endl;
+  else   geom->PrintGeometry();
+  
+  //Load Digits
+  rl->LoadSDigits("EMCAL");
+  
+  //Get maximum number of events
+  Int_t maxevent =  rl->GetNumberOfEvents();
+  cout<<"Number of events "<<maxevent<<endl;
+  //maxevent = 10 ;
+  
+  Int_t iEvent  = -1 ;
+  Float_t amp   = -1 ;
+  Float_t time  = -1 ;
+  Int_t id      = -1 ;
+  Int_t iSupMod =  0 ;
+  Int_t iTower  =  0 ;
+  Int_t iIphi   =  0 ;
+  Int_t iIeta   =  0 ;
+  Int_t iphi    =  0 ;
+  Int_t ieta    =  0 ;
+  
+  AliEMCALDigit * dig;
+  
+  for ( iEvent=0; iEvent<maxevent; iEvent++)
+    {
+      cout <<  " ======> Event " << iEvent << endl ;
+      //Load Event
+      rl->GetEvent(iEvent);
+      
+      //Fill array of digits
+      TClonesArray *digits = emcalLoader->SDigits();
+      
+      //Get digits from the list      
+      for(Int_t idig = 0; idig< digits->GetEntries();idig++){
+       //cout<<">> idig "<<idig<<endl;
+       dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
+       
+       if(dig != 0){
+         id   = dig->GetId() ; //cell (digit) label
+         amp  = dig->GetAmplitude(); //amplitude in cell (digit)
+         time = dig->GetTime();//time of creation of digit after collision
+         
+         cout<<"Cell ID "<<id<<" Amp "<<amp<<endl;//" time "<<time<<endl;
+         
+         //Geometry methods  
+         if(geom){
+           geom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta); 
+           //Gives SuperModule and Tower numbers
+           geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
+                                             iIphi, iIeta,iphi,ieta);
+           //Gives label of cell in eta-phi position per each supermodule
+           cout<< "SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
+               <<"; Phi "<<iIphi<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
+         }
+       }
+       else
+         cout<<"Digit pointer 0x0"<<endl;
+      }
+      
+    }
+
+
+}
+

diff --git a/EMCAL/macros/TestEMCALSimulation.C b/EMCAL/macros/TestEMCALSimulation.C
index 94ada2a7d87505bb1e8882f8c1bb22e4f80f23ed..d581b709c33f6a0cbc7a097cfce378c41ccc9be1 100644 (file)
--- a/EMCAL/macros/TestEMCALSimulation.C
+++ b/EMCAL/macros/TestEMCALSimulation.C
@@ -12,7 +12,11 @@ void TestEMCALSimulation(Int_t nev =10, Bool_t raw = kFALSE) {
   AliSimulation simulator;
   simulator.SetConfigFile("Config.C");
   simulator.SetMakeSDigits("EMCAL");
-  simulator.SetMakeDigits("EMCAL");
+  simulator.SetMakeDigits ("EMCAL");
+
+  //simulator.SetRunGeneration(kFALSE); // Generate or not particles
+  //simulator.SetRunSimulation(kFALSE); // Generate or not HITS (detector response) or not, start from SDigits
+
   if(raw)  simulator.SetWriteRawData("EMCAL","raw.root",kTRUE);
 
   //OCDB settings
@@ -20,13 +24,22 @@ void TestEMCALSimulation(Int_t nev =10, Bool_t raw = kFALSE) {
   simulator.SetSpecificStorage("GRP/GRP/Data",
                                Form("local://%s",gSystem->pwd()));
 
+  // In case of anchoring MC, comment previous OCDB lines
+  // select the appropriate year
+  //simulator.SetDefaultStorage("alien://Folder=/alice/data/2011/OCDB");
+  //simulator.UseVertexFromCDB();
+  //simulator.UseMagFieldFromGRP();
+
   //Avoid the HLT to run
   simulator.SetRunHLT("");
 
+  //Avoid QA
+  simulator.SetRunQA(":");
+  
   TStopwatch timer;
   timer.Start();
 
-//  simulator.SetRunNumber(140234);
+  simulator.SetRunNumber(159582); // LHC11d run
   
   simulator.Run(nev);
 

diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx
index aac7806b4d12a19293a5c33beaf0b5de53f80f75..e6639177642a62d682f8950e21bd61168fe519a9 100644 (file)
--- a/JETAN/AliAnalysisTaskJetCluster.cxx
+++ b/JETAN/AliAnalysisTaskJetCluster.cxx
@@ -106,6 +106,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
   fEventSelection(kFALSE),     
   fRequireVZEROAC(kFALSE),     
   fRequireTZEROvtx(kFALSE),
+  fUseHFcuts(kFALSE),
   fFilterMask(0),
   fFilterMaskBestPt(0),
   fFilterType(0),
@@ -275,6 +276,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fEventSelection(kFALSE),
   fRequireVZEROAC(kFALSE),     
   fRequireTZEROvtx(kFALSE), 
+  fUseHFcuts(kFALSE),
   fFilterMask(0),
   fFilterMaskBestPt(0),
   fFilterType(0),
@@ -1889,6 +1891,23 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
          if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
          continue;
        }
+
+       // heavy flavor jets
+       if(fFilterMask==528 && fUseHFcuts){
+          Double_t ntpcClus = tr->GetTPCNcls();
+          Double_t trPt=tr->Pt();
+         TFormula NTPCClsCut("f1NClustersTPCLinearPtDep","70.+30./20.*x");
+       
+         if (trPt <= 20. && (ntpcClus < NTPCClsCut.Eval(trPt))) continue;
+         else if (trPt > 20. && ntpcClus < 100) continue;
+
+         Int_t idtr1 = tr->GetID(); 
+         Bool_t flagSame = kFALSE;
+         AvoidDoubleCountingHF(aod, idtr1, flagSame);
+         if (flagSame) continue; 
+       }
+       // 
+
         if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
         if (fApplySharedClusterCut) {
            Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
@@ -2090,8 +2109,22 @@ Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *pa
        }                       
 return count;  
 }
+void AliAnalysisTaskJetCluster::AvoidDoubleCountingHF(AliAODEvent *aod, Int_t idtr1,  Bool_t &fFlagSameTr){
+  
+       Bool_t sametr = kFALSE;
 
+     for(int jt = 0;jt < aod->GetNumberOfTracks();++jt){
 
+       AliAODTrack *tr2 = aod->GetTrack(jt);
+       Int_t idtr2 = tr2->GetID();
+       if (idtr2>-1) continue;
+
+       idtr2=-1*idtr2;
+       if (idtr1==idtr2-1) sametr = kTRUE;
+       
+     }
+     fFlagSameTr = sametr;
+}
 
 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
 

diff --git a/JETAN/AliAnalysisTaskJetCluster.h b/JETAN/AliAnalysisTaskJetCluster.h
index 849f83708f1672acd4a751c3d4d69cedba732353..8424bd5b5ce068c107e1ce53650132bde1ed5a75 100644 (file)
--- a/JETAN/AliAnalysisTaskJetCluster.h
+++ b/JETAN/AliAnalysisTaskJetCluster.h
@@ -106,12 +106,12 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     virtual void SetFixedEfficiency(Double_t eff) {fEfficiencyFixed = eff;}
     virtual void SetRequireT0vtx(Bool_t b = true){fRequireTZEROvtx = b;}
     virtual void SetRequireV0AC(Bool_t b = true){fRequireVZEROAC = b;}
-
+    virtual void SetUseHFcuts(Bool_t b = true){fUseHFcuts = b;}
     Double_t GetMomentumSmearing(Int_t cat, Double_t pt);
     void FitMomentumResolution();
 
 
-    // for Fast Jet
+    // for Fast Jet 
     fastjet::JetAlgorithm        GetAlgorithm()         const {return fAlgorithm;}
     fastjet::Strategy            GetStrategy()          const {return fStrategy;}
     fastjet::RecombinationScheme GetRecombScheme()      const {return fRecombScheme;}
@@ -158,6 +158,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
 
     Int_t GetListOfTracks(TList *list,Int_t type);
        Int_t AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca);
+    void AvoidDoubleCountingHF(AliAODEvent *aod, Int_t idtr1, Bool_t &fFlagSameTr);
        
     AliAODEvent     *fAOD;                // ! where we take the jets from can be input or output AOD
     AliAODExtension *fAODExtension;       // ! AOD extension in case we write a non-sdt branch to a separate file and the aod is standard
@@ -168,6 +169,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     Bool_t        fEventSelection;        // use the event selection of this task, otherwise analyse all
     Bool_t        fRequireVZEROAC;        // switch to require V0 AC
     Bool_t        fRequireTZEROvtx;       // switch to require T0 vtx
+    Bool_t        fUseHFcuts;       // switch to require T0 vtx
     UInt_t        fFilterMask;            // filter bit for slecected tracks
     UInt_t        fFilterMaskBestPt;      // filter bit to mark jets with high quality leading tracks
 

diff --git a/PWGLF/STRANGENESS/Cascades/AliAnalysisTaskStrangenessVsMultiplicityMC.cxx b/PWGLF/STRANGENESS/Cascades/AliAnalysisTaskStrangenessVsMultiplicityMC.cxx
index 47b39941b03e22609c0d8cd1f4ba939f6fd0f6e5..ccb0096f94612ac6f035ee9bdacb04e12be64ced 100644 (file)
--- a/PWGLF/STRANGENESS/Cascades/AliAnalysisTaskStrangenessVsMultiplicityMC.cxx
+++ b/PWGLF/STRANGENESS/Cascades/AliAnalysisTaskStrangenessVsMultiplicityMC.cxx
@@ -609,7 +609,7 @@ void AliAnalysisTaskStrangenessVsMultiplicityMC::UserCreateOutputObjects()
   if(! fHistPtVsRefMultEta8_GenXiPlus ) {
     fHistPtVsRefMultEta8_GenXiPlus     = new TH2D( "fHistPtVsRefMultEta8_GenXiPlus",        "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenXiPlus);    }
   if(! fHistPtVsRefMultEta8_GenOmegaMinus ) {
-    fHistPtVsRefMultEta5_GenOmegaMinus    = new TH2D( "fHistPtVsRefMultEta8_GenOmegaMinus", "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenOmegaMinus);    }
+    fHistPtVsRefMultEta8_GenOmegaMinus    = new TH2D( "fHistPtVsRefMultEta8_GenOmegaMinus", "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenOmegaMinus);    }
   if(! fHistPtVsRefMultEta8_GenOmegaPlus ) {
     fHistPtVsRefMultEta8_GenOmegaPlus     = new TH2D( "fHistPtVsRefMultEta8_GenOmegaPlus",  "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenOmegaPlus);    }