Raw data decoder is migrated from AliCaloRawStream to AliCaloRawStreamV3
authorkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 3 Jul 2009 17:03:50 +0000 (17:03 +0000)
committerkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 3 Jul 2009 17:03:50 +0000 (17:03 +0000)
16 files changed:
PHOS/AliPHOS.cxx
PHOS/AliPHOSQADataMakerRec.cxx
PHOS/AliPHOSRawDigiProducer.cxx
PHOS/AliPHOSRawDigiProducer.h
PHOS/AliPHOSRawFitterv0.cxx [new file with mode: 0644]
PHOS/AliPHOSRawFitterv0.h [new file with mode: 0644]
PHOS/AliPHOSRawFitterv1.cxx [new file with mode: 0644]
PHOS/AliPHOSRawFitterv1.h [new file with mode: 0644]
PHOS/AliPHOSRawFitterv2.cxx [new file with mode: 0644]
PHOS/AliPHOSRawFitterv2.h [new file with mode: 0644]
PHOS/AliPHOSRecoParam.cxx
PHOS/AliPHOSRecoParam.h
PHOS/AliPHOSReconstructor.cxx
PHOS/CMake_libPHOSbase.txt
PHOS/PHOSbaseLinkDef.h
PHOS/libPHOSbase.pkg

index 7cbced5..efe9516 100644 (file)
@@ -94,8 +94,7 @@ class TFile;
 #include "AliPHOSCalibData.h"
 #include "AliPHOSPulseGenerator.h"
 #include "AliDAQ.h"
-#include "AliPHOSRawDecoder.h"
-#include "AliPHOSRawDecoderv1.h"
+#include "AliPHOSRawFitterv0.h"
 #include "AliPHOSCalibData.h"
 #include "AliPHOSRawDigiProducer.h"
 #include "AliPHOSQAChecker.h"
@@ -658,15 +657,15 @@ Bool_t AliPHOS::Raw2SDigits(AliRawReader* rawReader)
     mapping[i] = (AliAltroMapping*)maps->At(i);
   }
 
-  AliPHOSRawDecoderv1 dc(rawReader,mapping);    
+  AliPHOSRawFitterv0 fitter;
 
-  dc.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals());
-  dc.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset());
-  dc.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold());
+  fitter.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals());
+  fitter.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset());
+  fitter.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold());
 
-  AliPHOSRawDigiProducer pr;
+  AliPHOSRawDigiProducer pr(rawReader,mapping);
   pr.SetSampleQualityCut(AliPHOSSimParam::GetInstance()->GetEMCSampleQualityCut());     
-  pr.MakeDigits(sdigits,&dc);   
+  pr.MakeDigits(sdigits,&fitter);
                 
   Int_t bufferSize = 32000 ;    
   // TBranch * sdigitsBranch = tree->Branch("PHOS",&sdigits,bufferSize);        
index 934f8af..54c821e 100644 (file)
 // --- Standard library ---
 
 // --- AliRoot header files ---
+#include "AliCaloRawStreamV3.h"
 #include "AliESDCaloCluster.h"
 #include "AliESDEvent.h"
 #include "AliLog.h"
 #include "AliPHOSQADataMakerRec.h"
-#include "AliQAChecker.h"
 #include "AliPHOSDigit.h" 
 #include "AliPHOSCpvRecPoint.h" 
 #include "AliPHOSEmcRecPoint.h" 
+#include "AliPHOSRawFitterv0.h"
+#include "AliPHOSRawFitterv1.h"
+#include "AliPHOSRawFitterv2.h"
+#include "AliPHOSReconstructor.h"
 #include "AliPHOSRecParticle.h" 
 #include "AliPHOSTrackSegment.h" 
-#include "AliPHOSRawDecoder.h"
-#include "AliPHOSRawDecoderv1.h"
-#include "AliPHOSRawDecoderv2.h"
-#include "AliPHOSReconstructor.h"
+#include "AliQAChecker.h"
+#include "AliRawReader.h"
 
 ClassImp(AliPHOSQADataMakerRec)
            
@@ -374,111 +376,153 @@ void AliPHOSQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   //Fill prepared histograms with Raw digit properties
 
   rawReader->Reset() ;
-  AliPHOSRawDecoder * decoder ;
-  if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0)
-    decoder=new AliPHOSRawDecoderv1(rawReader);
-  else
-    if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v2")==0)
-      decoder=new AliPHOSRawDecoderv2(rawReader);
+
+  const TObjArray* maps = AliPHOSRecoParam::GetMappings();
+  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
+
+  AliAltroMapping *mapping[4];
+  for(Int_t i = 0; i < 4; i++) {
+    mapping[i] = (AliAltroMapping*)maps->At(i);
+  }
+
+  AliCaloRawStreamV3 *fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping);
+
+  AliPHOSRawFitterv0 * fitter ;
+  if     (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0)
+    fitter=new AliPHOSRawFitterv1();
+  else if(strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0)
+      fitter=new AliPHOSRawFitterv2();
     else
-      decoder=new AliPHOSRawDecoder(rawReader);
-  //decoder->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
-  decoder->SubtractPedestals(kTRUE);
+      fitter=new AliPHOSRawFitterv0();
+  fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
+//   fitter->SubtractPedestals(kTRUE);
   Double_t lgEtot=0. ;
   Double_t hgEtot=0. ;
-  Int_t lgNtot=0 ;
-  Int_t hgNtot=0 ;
-
-  while (decoder->NextDigit()) {
-   Int_t module  = decoder->GetModule() ;
-   Int_t row     = decoder->GetRow() ;
-   Int_t col     = decoder->GetColumn() ;
-   Double_t time = decoder->GetTime() ;
-   Double_t energy  = decoder->GetEnergy() ;
-   Bool_t lowGain = decoder->IsLowGain();
-   if (lowGain) {
-     //if(GetRecoParam()->EMCSubtractPedestals())
-       GetRawsData(kLGpedRMS)->Fill(decoder->GetPedestalRMS()) ;
-     if(energy<2.)
-       continue ;
-     switch(module){
-        case 1: GetRawsData(kLGmod1)->Fill(row-0.5,col-0.5) ; break ;
-        case 2: GetRawsData(kLGmod2)->Fill(row-0.5,col-0.5) ; break ;
-        case 3: GetRawsData(kLGmod3)->Fill(row-0.5,col-0.5) ; break ;
-        case 4: GetRawsData(kLGmod4)->Fill(row-0.5,col-0.5) ; break ;
-        case 5: GetRawsData(kLGmod5)->Fill(row-0.5,col-0.5) ; break ;
-     }
-     switch (module){
-        case 1: ((TH2F*)GetRawsData(kLGpedRMSMod1))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-        case 2: ((TH2F*)GetRawsData(kLGpedRMSMod2))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-        case 3: ((TH2F*)GetRawsData(kLGpedRMSMod3))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-        case 4: ((TH2F*)GetRawsData(kLGpedRMSMod4))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-        case 5: ((TH2F*)GetRawsData(kLGpedRMSMod5))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-     }
-     //if quality was evaluated, fill histo
-     if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0){
-       switch (module){
-         case 1: ((TH2F*)GetRawsData(kLGqualMod1))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 2: ((TH2F*)GetRawsData(kLGqualMod2))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 3: ((TH2F*)GetRawsData(kLGqualMod3))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 4: ((TH2F*)GetRawsData(kLGqualMod4))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 5: ((TH2F*)GetRawsData(kLGqualMod5))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-       }
-     }                                  
-     GetRawsData(kNmodLG)->Fill(module) ;
-     GetRawsData(kLGtime)->Fill(time) ; 
-     GetRawsData(kSpecLG)->Fill(energy) ;    
-     lgEtot+=energy ;
-     lgNtot++ ;   
-   } else {        
-     //if this isnon-ZS run - fill pedestal RMS
-     //if(GetRecoParam()->EMCSubtractPedestals())
-       GetRawsData(kHGpedRMS)->Fill(decoder->GetPedestalRMS()) ;
-     if(energy<8.)
-       continue ;
-     switch (module){
-         case 1: GetRawsData(kHGmod1)->Fill(row-0.5,col-0.5) ; break ;
-         case 2: GetRawsData(kHGmod2)->Fill(row-0.5,col-0.5) ; break ;
-         case 3: GetRawsData(kHGmod3)->Fill(row-0.5,col-0.5) ; break ;
-         case 4: GetRawsData(kHGmod4)->Fill(row-0.5,col-0.5) ; break ;
-         case 5: GetRawsData(kHGmod5)->Fill(row-0.5,col-0.5) ; break ;
-     }
-     switch (module){
-         case 1: ((TH2F*)GetRawsData(kHGpedRMSMod1))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-         case 2: ((TH2F*)GetRawsData(kHGpedRMSMod2))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-         case 3: ((TH2F*)GetRawsData(kHGpedRMSMod3))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-         case 4: ((TH2F*)GetRawsData(kHGpedRMSMod4))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-         case 5: ((TH2F*)GetRawsData(kHGpedRMSMod5))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ;
-     }               
-     //if quality was evaluated, fill histo
-     if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0){
-       switch (module){
-         case 1: ((TH2F*)GetRawsData(kHGqualMod1))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 2: ((TH2F*)GetRawsData(kHGqualMod2))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 3: ((TH2F*)GetRawsData(kHGqualMod3))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 4: ((TH2F*)GetRawsData(kHGqualMod4))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-         case 5: ((TH2F*)GetRawsData(kHGqualMod5))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ;
-       }
-     }
-     GetRawsData(kNmodHG)->Fill(module) ; 
-     GetRawsData(kHGtime)->Fill(time) ;  
-     GetRawsData(kSpecHG)->Fill(energy) ;
-     hgEtot+=energy ; 
-     hgNtot++ ;  
-   }                 
-  }                    
-  delete decoder;
+  Int_t    lgNtot=0 ;
+  Int_t    hgNtot=0 ;
+
+  while (fRawStream->NextDDL()) {   // !!!!!!!!!!! YK
+    while (fRawStream->NextChannel()) {
+      Int_t module   = fRawStream->GetModule();
+      Int_t cellX    = fRawStream->GetCellX();
+      Int_t cellZ    = fRawStream->GetCellZ();
+      Int_t caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU
+
+      Int_t nBunches = 0;
+      while (fRawStream->NextBunch()) {
+       nBunches++;
+       if (nBunches > 1) continue;
+       const UShort_t *sig = fRawStream->GetSignals();
+       Int_t sigStart      = fRawStream->GetStartTimeBin();
+       Int_t sigLength     = fRawStream->GetBunchLength();
+       fitter->SetSamples(sig,sigStart,sigLength);
+      } // End of NextBunch()
+
+      fitter->SetNBunches(nBunches);
+      fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
+      fitter->Eval();
+      
+      Double_t energy = fitter->GetEnergy() ; 
+      Double_t time   = fitter->GetTime() ;
+
+      if (caloFlag == 0) { // LG
+       //if(GetRecoParam()->EMCSubtractPedestals())
+       GetRawsData(kLGpedRMS)->Fill(fitter->GetPedestalRMS()) ;
+       if(energy<2.)
+         continue ;
+       switch(module){
+        case 1: GetRawsData(kLGmod1)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+        case 2: GetRawsData(kLGmod2)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+        case 3: GetRawsData(kLGmod3)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+        case 4: GetRawsData(kLGmod4)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+        case 5: GetRawsData(kLGmod5)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       }
+       switch (module){
+        case 1: ((TH2F*)GetRawsData(kLGpedRMSMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+        case 2: ((TH2F*)GetRawsData(kLGpedRMSMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+        case 3: ((TH2F*)GetRawsData(kLGpedRMSMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+        case 4: ((TH2F*)GetRawsData(kLGpedRMSMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+        case 5: ((TH2F*)GetRawsData(kLGpedRMSMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       }
+       //if quality was evaluated, fill histo
+       if(strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0){
+         switch (module){
+         case 1: ((TH2F*)GetRawsData(kLGqualMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 2: ((TH2F*)GetRawsData(kLGqualMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 3: ((TH2F*)GetRawsData(kLGqualMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 4: ((TH2F*)GetRawsData(kLGqualMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 5: ((TH2F*)GetRawsData(kLGqualMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         }
+       }                                  
+       GetRawsData(kNmodLG)->Fill(module) ;
+       GetRawsData(kLGtime)->Fill(time) ; 
+       GetRawsData(kSpecLG)->Fill(energy) ;    
+       lgEtot+=energy ;
+       lgNtot++ ;   
+      } else if (caloFlag == 1) { // HG        
+       //if this isnon-ZS run - fill pedestal RMS
+       //if(GetRecoParam()->EMCSubtractPedestals())
+       GetRawsData(kHGpedRMS)->Fill(fitter->GetPedestalRMS()) ;
+       if(energy<8.)
+         continue ;
+       switch (module){
+       case 1: GetRawsData(kHGmod1)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       case 2: GetRawsData(kHGmod2)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       case 3: GetRawsData(kHGmod3)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       case 4: GetRawsData(kHGmod4)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       case 5: GetRawsData(kHGmod5)->Fill(cellX-0.5,cellZ-0.5) ; break ;
+       }
+       switch (module){
+       case 1: ((TH2F*)GetRawsData(kHGpedRMSMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       case 2: ((TH2F*)GetRawsData(kHGpedRMSMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       case 3: ((TH2F*)GetRawsData(kHGpedRMSMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       case 4: ((TH2F*)GetRawsData(kHGpedRMSMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       case 5: ((TH2F*)GetRawsData(kHGpedRMSMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ;
+       }               
+       //if quality was evaluated, fill histo
+       if(strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0){
+         switch (module){
+         case 1: ((TH2F*)GetRawsData(kHGqualMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 2: ((TH2F*)GetRawsData(kHGqualMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 3: ((TH2F*)GetRawsData(kHGqualMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 4: ((TH2F*)GetRawsData(kHGqualMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         case 5: ((TH2F*)GetRawsData(kHGqualMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ;
+         }
+       }
+       GetRawsData(kNmodHG)->Fill(module) ; 
+       GetRawsData(kHGtime)->Fill(time) ;  
+       GetRawsData(kSpecHG)->Fill(energy) ;
+       hgEtot+=energy ; 
+       hgNtot++ ;  
+      }
+    }  // End of NextChannel
+  } // End of NextDDL
+  delete fitter;
+
   GetRawsData(kEtotLG)->Fill(lgEtot) ; 
-  TParameter<double> * p = dynamic_cast<TParameter<double>*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kEtotLG)->GetName()))) ; 
+  TParameter<double> * p;
+  p = dynamic_cast<TParameter<double>*>(GetParameterList()->
+                                       FindObject(Form("%s_%s_%s", GetName(), 
+                                                       AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), 
+                                                       GetRawsData(kEtotLG)->GetName()))) ; 
   if (p) p->SetVal(lgEtot) ; 
   GetRawsData(kEtotHG)->Fill(hgEtot) ;  
-  p = dynamic_cast<TParameter<double>*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kEtotHG)->GetName()))) ; 
+  p = dynamic_cast<TParameter<double>*>(GetParameterList()->
+                                       FindObject(Form("%s_%s_%s", GetName(), 
+                                                       AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), 
+                                                       GetRawsData(kEtotHG)->GetName()))) ; 
   if (p) p->SetVal(hgEtot) ; 
   GetRawsData(kNtotLG)->Fill(lgNtot) ;
-  p = dynamic_cast<TParameter<double>*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kNtotLG)->GetName()))) ; 
+  p = dynamic_cast<TParameter<double>*>(GetParameterList()->
+                                       FindObject(Form("%s_%s_%s", GetName(), 
+                                                       AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), 
+                                                       GetRawsData(kNtotLG)->GetName()))) ; 
   if (p) p->SetVal(lgNtot) ; 
   GetRawsData(kNtotHG)->Fill(hgNtot) ;
-  p = dynamic_cast<TParameter<double>*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kNtotHG)->GetName()))) ; 
+  p = dynamic_cast<TParameter<double>*>(GetParameterList()->
+                                       FindObject(Form("%s_%s_%s", 
+                                                       GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), 
+                                                       GetRawsData(kNtotHG)->GetName()))) ; 
   if (p) p->SetVal(hgNtot) ; 
 }
 
index ea39f00..5ee95aa 100644 (file)
@@ -16,7 +16,7 @@
 /* $Id$ */
 
 //This class produces PHOS digits of one event
-//using AliPHOSRawDecoder. 
+//using AliPHOSRawFitter. 
 //
 //   For example:
 //   TClonesArray *digits = new TClonesArray("AliPHOSDigit",100);
 
 // --- AliRoot header files ---
 #include "AliPHOSRawDigiProducer.h"
-#include "AliPHOSRawDecoder.h"
+#include "AliPHOSRawFitterv0.h"
 #include "AliPHOSGeometry.h"
 #include "AliPHOSDigit.h"
 #include "AliPHOSCalibData.h"
 #include "AliPHOSPulseGenerator.h"
+#include "AliCaloRawStreamV3.h"
 #include "AliLog.h"
 
 ClassImp(AliPHOSRawDigiProducer)
@@ -53,7 +54,9 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
   fSampleQualityCut(1.),
   fEmcCrystals(0),
   fGeom(0),
-  fPulseGenerator(0)
+  fPulseGenerator(0),
+  fRawReader(0),
+  fRawStream(0)
 {
   // Default constructor
 
@@ -65,7 +68,32 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
   GetCalibrationParameters() ; 
 
 }
-//--------------------------------------------------------------------------------------           
+//--------------------------------------------------------------------------------------
+AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader,
+                                              AliAltroMapping **mapping):
+  TObject(),
+  fEmcMinE(0.),
+  fCpvMinE(0.),
+  fSampleQualityCut(1.),
+  fEmcCrystals(0),
+  fGeom(0),
+  fPulseGenerator(0),
+  fRawReader(rawReader),
+  fRawStream(0)
+{
+  // Default constructor
+
+  fGeom=AliPHOSGeometry::GetInstance() ;
+  if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
+
+  fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
+  fPulseGenerator = new AliPHOSPulseGenerator();
+  GetCalibrationParameters() ; 
+
+  fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping);
+
+}
+//--------------------------------------------------------------------------------------
 AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
   TObject(),
   fEmcMinE(0.),
@@ -73,7 +101,9 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp)
   fSampleQualityCut(1.),
   fEmcCrystals(0),
   fGeom(0),
-  fPulseGenerator(0)
+  fPulseGenerator(0),
+  fRawReader(0),
+  fRawStream(0)
 {                                                          
   // Copy constructor
 
@@ -100,23 +130,24 @@ AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiP
   fPulseGenerator = new AliPHOSPulseGenerator();
   return  *this;
 } 
-//--------------------------------------------------------------------------------------                                                   
+//--------------------------------------------------------------------------------------
 AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
 {
-  // Desctructor
+  // Destructor
   if(fPulseGenerator) delete fPulseGenerator ;
   fPulseGenerator=0 ;
+  delete fRawStream;
 }
 //--------------------------------------------------------------------------------------
-void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder) 
+void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter) 
 {
   //Makes the job.
-  //TClonesArray *digits and raw data decoder should be provided by calling function.
+  //TClonesArray *digits and raw data fitter should be provided by calling function.
 
   digits->Clear();
  
-  Int_t    iDigit   = 0 ;
-  Int_t relId[4], absId =0;
+  Int_t iDigit=0 ;
+  Int_t relId[4], absId=-1, caloFlag=-1;
 
   const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch. 
   const Double_t highLowDiff=2.; //Maximal difference between High and Low channels in LG adc channels 
@@ -125,95 +156,109 @@ void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder*
   TClonesArray tmpLG("AliPHOSDigit",10000) ;
   Int_t ilgDigit=0 ;
 
-  //Let decoder subtract pedestals in case of ZS
-  decoder->SetCalibData(fgCalibData) ;
+  //Let fitter subtract pedestals in case of ZS
+  fitter->SetCalibData(fgCalibData) ;
   
-  while (decoder->NextDigit()) {
-
-    Double_t energy=decoder->GetEnergy() ; 
-    if(energy<=baseLine) //in ADC channels
-      continue ;
-
-    //remove digits with bad shape. Decoder should calculate quality so that 
-    //in default case quality [0,1], while larger values of quality mean somehow 
-    //corrupted samples, 999 means obviously corrupted sample.
-    //It is difficult to fit samples with overflow (even setting cut on overflow values)
-    //because too few points are left to fit. So we do not evaluate samples with overflow
-    if(decoder->GetSampleQuality() > fSampleQualityCut && !(decoder->IsOverflow()))
-       continue ;
-
-    Bool_t lowGainFlag = decoder->IsLowGain();
-
-    relId[0] = decoder->GetModule();
-    relId[1] = 0;
-    relId[2] = decoder->GetRow();
-    relId[3] = decoder->GetColumn();
-    fGeom->RelToAbsNumbering(relId, absId);
-
-    Double_t time = decoder->GetTime() ;
-    time = CalibrateT(time,relId,lowGainFlag) ;
-
-    energy = CalibrateE(energy,relId,lowGainFlag) ;
-
-    if(energy <= 0.) 
-       continue;
-
-    if(lowGainFlag){
-      new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
-      ilgDigit++ ; 
-    }
-    else{ 
-      if(decoder->IsOverflow()) //Keep this digit to replace it by Low Gain later.
-                                //If there is no LogGain it wil be removed by cut on Min E
-        new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
-      else
-        new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
-      iDigit++;
-    }
-  }
-
-  //Now scan created LG and HG digits and keep only those which appeared in both lists 
-  //replace energy of HighGain digits only if there is overflow
-  //negative energy (overflow)
-  digits->Sort() ;
-  tmpLG.Sort() ;
-  Int_t iLG = 0;
-  Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
-
-  for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) { 
-    AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
-    if (!digHG) continue;
-    AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
-    while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
-      iLG++ ;
-      digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
-    }
-    absId=digHG->GetId() ;                                                                                                         
-    fGeom->AbsToRelNumbering(absId,relId) ;                                                                                              
+  while (fRawStream->NextDDL()) {
+    while (fRawStream->NextChannel()) {
+      relId[0] = fRawStream->GetModule();
+      relId[1] = 0;
+      relId[2] = fRawStream->GetCellX();
+      relId[3] = fRawStream->GetCellZ();
+      caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU
+      fGeom->RelToAbsNumbering(relId, absId);
+
+      Int_t nBunches = 0;
+      while (fRawStream->NextBunch()) {
+       nBunches++;
+       if (nBunches > 1) continue;
+       const UShort_t *sig = fRawStream->GetSignals();
+       Int_t sigStart  = fRawStream->GetStartTimeBin();
+       Int_t sigLength = fRawStream->GetBunchLength();
+       fitter->SetSamples(sig,sigStart,sigLength);
+      } // End of NextBunch()
+
+      fitter->SetNBunches(nBunches);
+      fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
+      fitter->Eval();
+      
+      Double_t energy = fitter->GetEnergy() ; 
+      Double_t time   = fitter->GetTime() ;
+      if(energy<=baseLine) //in ADC channels
+       continue ;
+
+      //remove digits with bad shape. Fitter should calculate quality so that 
+      //in default case quality [0,1], while larger values of quality mean somehow 
+      //corrupted samples, 999 means obviously corrupted sample.
+      //It is difficult to fit samples with overflow (even setting cut on overflow values)
+      //because too few points are left to fit. So we do not evaluate samples with overflow
+
+      if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
+       continue ;
+      
+//       energy = CalibrateE(energy,relId,lowGainFlag) ;
+//       time   = CalibrateT(time,relId,lowGainFlag) ;
+      
+      if(energy <= 0.) 
+       continue;
+      
+      if (caloFlag == AliCaloRawStreamV3::kLowGain) {
+       new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
+       ilgDigit++ ; 
+      }
+      else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
+       if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
+         //If there is no LogGain it wil be removed by cut on Min E
+         new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
+       else
+         new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
+       iDigit++;
+      }
+    } // End of NextChannel()
+
+    //Now scan created LG and HG digits and keep only those which appeared in both lists 
+    //replace energy of HighGain digits only if there is overflow
+    //negative energy (overflow)
+    digits->Sort() ;
+    tmpLG.Sort() ;
+    Int_t iLG = 0;
+    Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
+    
+    for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) { 
+      AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
+      if (!digHG) continue;
+      AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
+      while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
+       iLG++ ;
+       digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
+      }
+      absId=digHG->GetId() ;
+      fGeom->AbsToRelNumbering(absId,relId) ;
  
-    if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
-      if(digHG->GetEnergy()<0.){ //This is overflow in HG
-        digHG->SetTime(digLG->GetTime()) ;
-        digHG->SetEnergy(digLG->GetEnergy()) ;
-      } 
-      else{ //Make approximate comparison of HG and LG energies
-        Double_t de = (digHG->GetEnergy()-digLG->GetEnergy()) ; 
-        if(TMath::Abs(de)>CalibrateE(double(highLowDiff),relId,1)){ //too strong difference, remove digit
-          digits->RemoveAt(iDig) ;
-        }
+      if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
+       if(digHG->GetEnergy()<0.){ //This is overflow in HG
+         digHG->SetTime(digLG->GetTime()) ;
+         digHG->SetEnergy(digLG->GetEnergy()) ;
+       }
+       else{ //Make approximate comparison of HG and LG energies
+         Double_t de = (digHG->GetEnergy()-digLG->GetEnergy()) ; 
+         if(TMath::Abs(de)>CalibrateE(double(highLowDiff),relId,1)){ //too strong difference, remove digit
+           digits->RemoveAt(iDig) ;
+         }
+       }
+      }
+      else{ //no pair - remove
+       // temporary fix for dead LG channels
+       if(relId[2]%2==1 && relId[3]%16==4) 
+         continue ;
+       if(digHG->GetEnergy()>CalibrateE(double(5.),relId,1)) //One can not always find LG with Amp<5 ADC ch.
+         digits->RemoveAt(iDig) ;                                                                                                            
       }
     }
-    else{ //no pair - remove
-      // temporary fix for dead LG channels
-      if(relId[2]%2==1 && relId[3]%16==4) 
-        continue ;
-      if(digHG->GetEnergy()>CalibrateE(double(5.),relId,1)) //One can not always find LG with Amp<5 ADC ch.
-        digits->RemoveAt(iDig) ;                                                                                                            
-    }
-  }
+  } // End of NextDDL()
 
   CleanDigits(digits) ;
-
+  
 }
 //____________________________________________________________________________
 Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
index 183c4db..352641c 100644 (file)
@@ -6,13 +6,17 @@
 /* $Id$ */
 
 //This class produces PHOS digits of one event
-//using AliPHOSRawDecoder. See cxx source for use case.
+//using AliPHOSRawFitter. See cxx source for use case.
 
-class AliPHOSRawDecoder;
 class AliPHOSCalibData ;
 class AliPHOSDigit ;
 class AliPHOSGeometry ;
 class AliPHOSPulseGenerator;
+class AliRawReader;
+class AliCaloRawStreamV3;
+class AliPHOSRawFitterv0;
+
+#include "AliAltroMapping.h"
 #include "TObject.h"
 
 class AliPHOSRawDigiProducer: public TObject {
@@ -20,12 +24,13 @@ class AliPHOSRawDigiProducer: public TObject {
 public:
 
   AliPHOSRawDigiProducer() ;
+  AliPHOSRawDigiProducer(AliRawReader *rawReader, AliAltroMapping **mapping = NULL);
   AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp);
   AliPHOSRawDigiProducer& operator= (const AliPHOSRawDigiProducer &dp);
  
   virtual ~AliPHOSRawDigiProducer(); 
 
-  void MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder);
+  void MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter);
 
   void SetEmcMinAmp(Float_t emcMin) { fEmcMinE=emcMin; }
   void SetCpvMinAmp(Float_t cpvMin) { fCpvMinE=cpvMin; }
@@ -48,10 +53,12 @@ private:
   Float_t fSampleQualityCut;         // Cut on sample shapes: 0: no samples; 1: default parameterization; 999: accept even obviously bad
   Int_t fEmcCrystals ;               //  number of EMC crystals
   AliPHOSGeometry * fGeom ;          //! PHOS geometry
-  static AliPHOSCalibData * fgCalibData ;   //! Calibration database if avalable
-  AliPHOSPulseGenerator * fPulseGenerator ; //! Class with pulse shape parameters
+  static AliPHOSCalibData * fgCalibData ;     //! Calibration database if avalable
+  AliPHOSPulseGenerator   * fPulseGenerator ; //! Class with pulse shape parameters
+  AliRawReader            * fRawReader;       //! Raw data reader
+  AliCaloRawStreamV3      * fRawStream;       //! Calorimeter decoder of ALTRO format
 
-  ClassDef(AliPHOSRawDigiProducer,5)
+  ClassDef(AliPHOSRawDigiProducer,6)
 };
 
 #endif
diff --git a/PHOS/AliPHOSRawFitterv0.cxx b/PHOS/AliPHOSRawFitterv0.cxx
new file mode 100644 (file)
index 0000000..e1d1125
--- /dev/null
@@ -0,0 +1,220 @@
+/**************************************************************************
+ * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+// This class extracts the signal parameters (energy, time, quality)
+// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
+// A coarse algorithm takes the energy as the maximum
+// sample, time is the first sample index, and the quality is the number
+// of bunches in the signal.
+// 
+//     AliPHOSRawFitterv0 *fitterv0=new AliPHOSRawFitterv0();
+//     fitterv0->SetSamples(sig,sigStart,sigLength);
+//     fitterv0->SetNBunches(nBunches);
+//     fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag);
+//     fitterv0->SetCalibData(fgCalibData) ;
+//     fitterv0->Eval();
+//     Double_t amplitude = fitterv0.GetEnergy();
+//     Double_t time      = fitterv0.GetTime();
+//     Bool_t   isLowGain = fitterv0.GetCaloFlag()==0;
+
+// Author: Yuri Kharlov
+
+// --- ROOT system ---
+#include "TArrayI.h"
+#include "TMath.h"
+#include "TObject.h"
+
+// --- AliRoot header files ---
+#include "AliPHOSRawFitterv0.h"
+#include "AliPHOSCalibData.h"
+#include "AliLog.h"
+
+ClassImp(AliPHOSRawFitterv0)
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv0::AliPHOSRawFitterv0():
+  TObject(),
+  fSignal(0),
+  fModule(0),
+  fCellX(0),
+  fCellZ(0),
+  fCaloFlag(0),
+  fStart(0),
+  fLength(0),
+  fNBunches(0),
+  fPedSubtract(kFALSE),
+  fEnergy(-111),
+  fTime(-111),
+  fQuality(0.),
+  fPedestalRMS(0.),
+  fAmpOffset(0),
+  fAmpThreshold(0),
+  fOverflow(kFALSE),
+  fCalibData(0)
+{
+  //Default constructor
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv0::~AliPHOSRawFitterv0()
+{
+  //Destructor
+  delete [] fSignal;
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ):
+  TObject(),
+  fSignal      (phosFitter.fSignal),
+  fModule      (phosFitter.fModule),
+  fCellX       (phosFitter.fCellX),
+  fCellZ       (phosFitter.fCellZ),
+  fCaloFlag    (phosFitter.fCaloFlag),
+  fStart       (phosFitter.fStart),
+  fLength      (phosFitter.fLength),
+  fNBunches    (phosFitter.fNBunches),
+  fPedSubtract (phosFitter.fPedSubtract),
+  fEnergy      (phosFitter.fEnergy),
+  fTime        (phosFitter.fTime),
+  fQuality     (phosFitter.fQuality),
+  fPedestalRMS (phosFitter.fPedestalRMS),
+  fAmpOffset   (phosFitter.fAmpOffset),
+  fAmpThreshold(phosFitter.fAmpThreshold),
+  fOverflow    (phosFitter.fOverflow),
+  fCalibData   (phosFitter.fCalibData)
+{
+  //Copy constructor
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter)
+{
+  //Assignment operator.
+
+  if(this != &phosFitter) {
+    fSignal       = phosFitter.fSignal;
+    fModule       = phosFitter.fModule;
+    fCellX        = phosFitter.fCellX;
+    fCellZ        = phosFitter.fCellZ;
+    fCaloFlag     = phosFitter.fCaloFlag;
+    fStart        = phosFitter.fStart;
+    fLength       = phosFitter.fLength;
+    fNBunches     = phosFitter.fNBunches;
+    fPedSubtract  = phosFitter.fPedSubtract;
+    fEnergy       = phosFitter.fEnergy;
+    fTime         = phosFitter.fTime;
+    fQuality      = phosFitter.fQuality;
+    fPedestalRMS  = phosFitter.fPedestalRMS;
+    fAmpOffset    = phosFitter.fAmpOffset;
+    fAmpThreshold = phosFitter.fAmpThreshold;
+    fOverflow     = phosFitter.fOverflow;
+    fCalibData    = phosFitter.fCalibData;
+  }
+
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+
+void AliPHOSRawFitterv0::SetSamples(const UShort_t *sig, Int_t sigStart, Int_t sigLength)
+{
+  // Set the sample array, its start and length in time bin units
+
+  fStart   = sigStart;
+  fLength  = sigLength;
+  fSignal  = new UShort_t[fLength];
+  for (Int_t i=0; i<fLength; i++) {
+    fSignal[i] = sig[i];
+  }
+}
+//-----------------------------------------------------------------------------
+
+void AliPHOSRawFitterv0::SetChannelGeo(const Int_t module, const Int_t cellX,
+                                    const Int_t cellZ,  const Int_t caloFlag)
+{
+  // Set geometry address of the channel
+  // (for a case if fitting parameters are different for different channels)
+  
+  fModule   = module;
+  fCellX    = cellX;
+  fCellZ    = cellZ;
+  fCaloFlag = caloFlag;
+}
+//-----------------------------------------------------------------------------
+
+Bool_t AliPHOSRawFitterv0::Eval()
+{
+  // Calculate signal parameters (energy, time, quality) from array of samples
+  // Energy is a maximum sample minus pedestal 9
+  // Time is the first time bin
+  // Signal overflows is there are at least 3 samples of the same amplitude above 900
+
+  fEnergy  = 0;
+  if (fNBunches > 1) {
+    fQuality = 1000;
+    return kTRUE;
+  }
+  
+  const Float_t kBaseLine   = 1.0;
+  const Int_t   kPreSamples = 10;
+
+  Float_t  pedMean   = 0;
+  Float_t  pedRMS    = 0;
+  Int_t    nPed      = 0;
+  UShort_t maxSample = 0;
+  Int_t    nMax      = 0;
+
+  for (Int_t i=0; i<fLength; i++) {
+    if (i<kPreSamples) {
+      nPed++;
+      pedMean += fSignal[i];
+      pedRMS  += fSignal[i]*fSignal[i] ;
+    }
+    if(fSignal[i] > maxSample) maxSample = fSignal[i];
+    if(fSignal[i] == maxSample) nMax++;
+  }
+  fEnergy = (Double_t)maxSample;
+  if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
+
+  if (fPedSubtract) {
+    if (nPed > 0) {
+      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
+      if(fPedestalRMS > 0.) 
+       fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
+      fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
+    }
+    else
+      return kFALSE;
+  }
+  else {
+    //take pedestals from DB
+    Double_t pedestal = (Double_t) fAmpOffset ;
+    if (fCalibData) {
+      Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
+      Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
+      pedestal += truePed - altroSettings ;
+    }
+    else{
+      AliWarning(Form("Can not read data from OCDB")) ;
+    }
+    fEnergy-=pedestal ;
+  }
+  if (fEnergy < kBaseLine) fEnergy = 0;
+  
+  return kTRUE;
+
+}
diff --git a/PHOS/AliPHOSRawFitterv0.h b/PHOS/AliPHOSRawFitterv0.h
new file mode 100644 (file)
index 0000000..631d1e7
--- /dev/null
@@ -0,0 +1,69 @@
+#ifndef ALIPHOSRAWFITTERV0_H
+#define ALIPHOSRAWFITTERV0_H
+/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                          */
+
+/* $Id: $ */
+
+// This class extracts the signal parameters (energy, time, quality)
+// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
+// A coarse algorithm is applied.
+
+class TArrayI;
+class AliPHOSCalibData;
+
+class AliPHOSRawFitterv0 : public TObject 
+{
+
+public:
+
+  AliPHOSRawFitterv0();
+  AliPHOSRawFitterv0(const AliPHOSRawFitterv0& rawFitterv0);
+  AliPHOSRawFitterv0& operator = (const AliPHOSRawFitterv0& rawFitterv0);
+  virtual ~AliPHOSRawFitterv0();
+
+  void SubtractPedestals(Bool_t subtract) {fPedSubtract  = subtract;}
+  void SetAmpOffset     (Int_t extPed=5)  {fAmpOffset    = extPed ;}
+  void SetAmpThreshold  (Int_t thr=5)     {fAmpThreshold = thr ;}
+  void SetSamples(const UShort_t *sig, Int_t sigStart, Int_t sigLength);
+  void SetNBunches(const Int_t nBunches) { fNBunches = nBunches; }
+  void SetChannelGeo(const Int_t module, const Int_t cellX,
+                    const Int_t cellZ,  const Int_t caloFlag);
+
+  virtual Bool_t Eval();
+  Double_t GetEnergy()        const { return fEnergy;      }
+  Double_t GetTime()          const { return fTime;        }
+  Double_t GetSignalQuality() const { return fQuality;     }
+  Double_t GetPedestalRMS()   const { return fPedestalRMS; }
+  Int_t    GetModule()        const { return fModule;      }
+  Int_t    GetCellX()         const { return fCellX;       }
+  Int_t    GetCellZ()         const { return fCellZ;       }
+  Int_t    GetCaloFlag()      const { return fCaloFlag;    }
+  Bool_t   IsOverflow()       const { return fOverflow;    }
+
+  void SetCalibData(AliPHOSCalibData * cdata){ fCalibData=cdata ;}
+
+protected:   
+  
+  UShort_t *fSignal;        // array of samples
+  Int_t    fModule;         // PHOS module number
+  Int_t    fCellX;          // cell number along X-axis
+  Int_t    fCellZ;          // cell number along Z-axis
+  Int_t    fCaloFlag;       // 0=LG, 1=HG, 2=TRU
+  Int_t    fStart;          // time bin of start signal
+  Int_t    fLength;         // signal length in time bins
+  Int_t    fNBunches;       // number of bunches in a signal
+  Bool_t   fPedSubtract;    // pedestals subtraction (kTRUE="yes")
+  Double_t fEnergy;         // "digit" energy
+  Double_t fTime;           // "digit" time
+  Double_t fQuality ;       // sample quality
+  Double_t fPedestalRMS;    // calciulated RMS of pedestal (non-ZS runs)
+  Int_t    fAmpOffset ;     // pedestal offset from ALTRO chips
+  Int_t    fAmpThreshold ;  // zero suppression threshold from ALTRO chips
+  Bool_t   fOverflow ;      // kTRUE is the signal overflows
+  AliPHOSCalibData * fCalibData ;   //! Calibration database if avalable
+
+  ClassDef(AliPHOSRawFitterv0,1)
+};
+
+#endif
diff --git a/PHOS/AliPHOSRawFitterv1.cxx b/PHOS/AliPHOSRawFitterv1.cxx
new file mode 100644 (file)
index 0000000..7359f79
--- /dev/null
@@ -0,0 +1,445 @@
+/**************************************************************************
+ * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+// This class extracts the signal parameters (energy, time, quality)
+// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
+// A fitting algorithm evaluates the energy and the time from Minuit minimization
+// 
+// Typical use case:
+//     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
+//     fitter->SetSamples(sig,sigStart,sigLength);
+//     fitter->SetNBunches(nBunches);
+//     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
+//     fitter->SetCalibData(fgCalibData) ;
+//     fitter->Eval();
+//     Double_t amplitude = fitter.GetEnergy();
+//     Double_t time      = fitter.GetTime();
+//     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
+
+// Author: Dmitri Peressounko (Oct.2008)
+// Modified: Yuri Kharlov (Jul.2009)
+
+// --- ROOT system ---
+#include "TArrayD.h"
+#include "TList.h"
+#include "TMath.h"
+#include "TMinuit.h"
+#include "TCanvas.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TF1.h"
+#include "TROOT.h"
+
+// --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSRawFitterv1.h"
+#include "AliPHOSPulseGenerator.h"
+
+ClassImp(AliPHOSRawFitterv1)
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv1::AliPHOSRawFitterv1():
+  AliPHOSRawFitterv0(),
+  fSampleParamsLow(0x0),
+  fSampleParamsHigh(0x0),
+  fToFit(0x0)
+{
+  //Default constructor.
+  if(!gMinuit) 
+    gMinuit = new TMinuit(100);
+  fSampleParamsHigh =new TArrayD(7) ;
+  fSampleParamsHigh->AddAt(2.174,0) ;
+  fSampleParamsHigh->AddAt(0.106,1) ;
+  fSampleParamsHigh->AddAt(0.173,2) ;
+  fSampleParamsHigh->AddAt(0.06106,3) ;
+  //last two parameters are pedestal and overflow
+  fSampleParamsLow=new TArrayD(7) ;
+  fSampleParamsLow->AddAt(2.456,0) ;
+  fSampleParamsLow->AddAt(0.137,1) ;
+  fSampleParamsLow->AddAt(2.276,2) ;
+  fSampleParamsLow->AddAt(0.08246,3) ;
+  fToFit = new TList() ;
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv1::~AliPHOSRawFitterv1()
+{
+  //Destructor.
+  if(fSampleParamsLow){
+    delete fSampleParamsLow ; 
+    fSampleParamsLow=0 ;
+  }
+  if(fSampleParamsHigh){
+    delete fSampleParamsHigh ;
+    fSampleParamsHigh=0;
+  }
+  if(fToFit){
+    delete fToFit ;
+    fToFit=0 ;
+  }
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv1::AliPHOSRawFitterv1(const AliPHOSRawFitterv1 &phosFitter ):
+  AliPHOSRawFitterv0(phosFitter), 
+  fSampleParamsLow(0x0),
+  fSampleParamsHigh(0x0),
+  fToFit(0x0)
+{
+  //Copy constructor.
+  fToFit = new TList() ;
+  fSampleParamsLow =new TArrayD(*(phosFitter.fSampleParamsLow)) ;
+  fSampleParamsHigh=new TArrayD(*(phosFitter.fSampleParamsHigh)) ;
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv1& AliPHOSRawFitterv1::operator = (const AliPHOSRawFitterv1 &phosFitter)
+{
+  //Assignment operator.
+
+  fToFit = new TList() ;
+  if(fSampleParamsLow){
+    fSampleParamsLow = phosFitter.fSampleParamsLow ;
+    fSampleParamsHigh= phosFitter.fSampleParamsHigh ;
+  }
+  else{
+    fSampleParamsLow =new TArrayD(*(phosFitter.fSampleParamsLow)) ; 
+    fSampleParamsHigh=new TArrayD(*(phosFitter.fSampleParamsHigh)) ;
+  }
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliPHOSRawFitterv1::Eval()
+{
+  //Extract an energy deposited in the crystal,
+  //crystal' position (module,column,row),
+  //time and gain (high or low).
+  //First collects sample, then evaluates it and if it has
+  //reasonable shape, fits it with Gamma2 function and extracts 
+  //energy and time.
+
+  if (fNBunches > 1) {
+    fQuality = 1000;
+    return kTRUE;
+  }
+
+  Float_t pedMean = 0;
+  Float_t pedRMS  = 0;
+  Int_t   nPed    = 0;
+  const Float_t kBaseLine   = 1.0;
+  const Int_t   kPreSamples = 10;
+  
+  TArrayI *fSamples = new TArrayI(fLength); // array of samples
+  TArrayI *fTimes   = new TArrayI(fLength); // array of times corresponding to samples
+  for (Int_t i=0; i<fLength; i++) {
+    if (i<kPreSamples) {
+      nPed++;
+      pedMean += fSignal[i];
+      pedRMS  += fSignal[i]*fSignal[i] ;
+    }
+    fSamples->AddAt(fSignal[i],i);
+    fTimes  ->AddAt(        i ,i);
+  }
+
+  Int_t    iBin     = fSamples->GetSize() ;
+  fEnergy = -111;
+  fQuality= 999. ;
+  const Float_t sampleMaxHG=102.332 ;  //maximal height of HG sample with given parameterization
+  const Float_t sampleMaxLG=277.196 ;  //maximal height of LG sample with given parameterization
+  const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
+  Double_t pedestal = 0;
+
+  if (fPedSubtract) {
+    if (nPed > 0) {
+      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
+      if(fPedestalRMS > 0.) 
+       fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
+      fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
+    }
+    else
+      return kFALSE;
+  }
+  else {
+    //take pedestals from DB
+    pedestal = (Double_t) fAmpOffset ;
+    if (fCalibData) {
+      Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
+      Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
+      pedestal += truePed - altroSettings ;
+    }
+    else{
+      AliWarning(Form("Can not read data from OCDB")) ;
+    }
+    fEnergy-=pedestal ;
+  }
+
+  if (fEnergy < kBaseLine) fEnergy = 0;
+  
+  //calculate time and energy
+  Int_t maxBin=0 ;
+  Int_t maxAmp=0 ;
+  Double_t aMean=0. ;
+  Double_t aRMS=0. ;
+  Double_t wts=0 ;
+  Int_t tStart = 0 ;
+  for(Int_t i=0; i<fSamples->GetSize(); i++){
+    if(fSamples->At(i) > pedestal){
+      Double_t de=fSamples->At(i)-pedestal ;
+      if(de>1.){
+       aMean+=de*i ;
+       aRMS+=de*i*i ;
+       wts+=de; 
+      }
+      if(de>2 && tStart==0) 
+       tStart=i ;
+      if(maxAmp<fSamples->At(i)){
+       maxBin=i ;
+       maxAmp=fSamples->At(i) ;
+      }
+    }
+  }
+  if(maxBin==fSamples->GetSize()-1){//bad "rising" sample
+    fEnergy=0. ;
+    fTime=-999.;
+    fQuality= 999. ;
+    return kTRUE ;
+  }
+  fEnergy=Double_t(maxAmp)-pedestal ;
+  fOverflow =0 ;  //look for plato on the top of sample
+  if(fEnergy>500 &&  //this is not fluctuation of soft sample
+     maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato
+    fOverflow = kTRUE ;
+  }
+      
+  if(wts>0){
+    aMean/=wts; 
+    aRMS=aRMS/wts-aMean*aMean;
+  }
+
+  //do not take too small energies
+  if(fEnergy < kBaseLine) 
+    fEnergy = 0;
+  
+  //do not test quality of too soft samples
+  if(fEnergy<maxEtoFit){
+    fTime=fTimes->At(tStart);
+    if(aRMS<2.) //sigle peak
+      fQuality=999. ;
+    else
+      fQuality= 0. ;
+    return kTRUE ;
+  }
+      
+  //IF sample has reasonable mean and RMS, try to fit it with gamma2
+  
+  gMinuit->mncler();                     // Reset Minuit's list of paramters
+  gMinuit->SetPrintLevel(-1) ;           // No Printout
+  gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ;  
+  // To set the address of the minimization function 
+  
+  fToFit->Clear("nodelete") ;
+  Double_t b=0,bmin=0,bmax=0 ;
+  if(fCaloFlag == 0){ // Low gain
+    fSampleParamsLow->AddAt(pedestal,4) ;
+    if(fOverflow)
+      fSampleParamsLow->AddAt(double(maxAmp),5) ;
+    else
+      fSampleParamsLow->AddAt(double(1023),5) ;
+    fSampleParamsLow->AddAt(double(iBin),6) ;
+    fToFit->AddFirst((TObject*)fSampleParamsLow) ; 
+    b=fSampleParamsLow->At(2) ;
+    bmin=0.5 ;
+    bmax=10. ;
+  }
+  else if(fCaloFlag == 1){ // High gain
+    fSampleParamsHigh->AddAt(pedestal,4) ;
+    if(fOverflow)
+      fSampleParamsHigh->AddAt(double(maxAmp),5) ;
+    else
+      fSampleParamsHigh->AddAt(double(1023),5);
+    fSampleParamsHigh->AddAt(double(iBin),6);
+    fToFit->AddFirst((TObject*)fSampleParamsHigh) ; 
+    b=fSampleParamsHigh->At(2) ;
+    bmin=0.05 ;
+    bmax=0.4 ;
+  }
+  fToFit->AddLast((TObject*)fSamples) ;
+  fToFit->AddLast((TObject*)fTimes) ;
+  
+  gMinuit->SetObjectFit((TObject*)fToFit) ;         // To tranfer pointer to UnfoldingChiSquare
+  Int_t ierflg ;
+  gMinuit->mnparm(0, "t0",  1.*tStart, 0.01, -500., 500., ierflg) ;
+  if(ierflg != 0){
+    //   AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
+    fEnergy=0. ;
+    fTime=-999. ;
+    fQuality=999 ;
+    return kTRUE ; //will scan further
+  }
+  Double_t amp0=0; 
+  if(fCaloFlag == 0) // Low gain
+    amp0=fEnergy/sampleMaxLG;
+  else if(fCaloFlag == 1) // High gain
+    amp0=fEnergy/sampleMaxHG;
+  
+  gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
+  if(ierflg != 0){
+    //   AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
+    fEnergy=0. ;
+    fTime=-999. ;
+    fQuality=999 ;
+    return kTRUE ; //will scan further
+  }
+  
+  gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ;
+  if(ierflg != 0){                                         
+    //        AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;  
+    fEnergy=0. ;           
+    fTime=-999. ;         
+    fQuality=999 ;       
+    return kTRUE ; //will scan further  
+  }             
+  
+  
+  Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
+  //  depends on it. 
+  Double_t p1 = 1.0 ;
+  Double_t p2 = 0.0 ;
+  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
+  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
+  //   gMinuit->SetMaxIterations(100);
+  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
+  
+  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
+  
+  Double_t err,t0err ;
+  Double_t t0,efit ;
+  gMinuit->GetParameter(0,t0, t0err) ;    
+  gMinuit->GetParameter(1,efit, err) ;    
+  
+  Double_t bfit, berr ;
+  gMinuit->GetParameter(2,bfit,berr) ;
+  
+  //Calculate total energy
+  //this isparameterization of depetendence of pulse height on parameter b
+  if(fCaloFlag == 0) // Low gain
+    efit*=99.54910 + 78.65038*bfit ;
+  else if(fCaloFlag == 1) // High gain
+    efit*=80.33109+128.6433*bfit ;
+  
+  if(efit<0. || efit > 10000.){
+    //set energy to previously found max
+    fTime=-999.;                                                                
+    fQuality=999 ;                                                              
+    return kTRUE;
+  }                                                                             
+  
+  //evaluate fit quality
+  Double_t fmin,fedm,errdef ;
+  Int_t npari,nparx,istat;
+  gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
+  fQuality=fmin/(fSamples->GetSize()-iBin) ;
+  //compare quality with some parameterization
+  if(fCaloFlag == 0) // Low gain
+    fQuality/=2.+0.002*fEnergy ;
+  else if(fCaloFlag == 1) // High gain
+    fQuality/=0.75+0.0025*fEnergy ;
+  
+  fEnergy=efit ;
+  fTime=t0-4.024*bfit ; //-10.402*bfit+4.669*bfit*bfit ; //Correction for 70 samples
+  
+  delete fSamples ;
+  delete fTimes ;
+  return kTRUE;
+}
+//_____________________________________________________________________________
+void AliPHOSRawFitterv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
+{
+  // Number of parameters, Gradient, Chi squared, parameters, what to do
+
+  TList * toFit= (TList*)gMinuit->GetObjectFit() ;
+  TArrayD * params=(TArrayD*)toFit->At(0) ; 
+  TArrayI * samples = (TArrayI*)toFit->At(1) ;
+  TArrayI * times = (TArrayI*)toFit->At(2) ;
+
+  fret = 0. ;     
+  if(iflag == 2)
+    for(Int_t iparam = 0 ; iparam < 3 ; iparam++)    
+      Grad[iparam] = 0 ; // Will evaluate gradient
+  
+  Double_t t0=x[0] ;
+  Double_t en=x[1] ;
+  Double_t b=x[2] ;
+  Double_t n=params->At(0) ;
+  Double_t alpha=params->At(1) ;
+  Double_t beta=params->At(3) ;
+  Double_t ped=params->At(4) ;
+
+  Double_t overflow=params->At(5) ;
+  Int_t iBin = (Int_t) params->At(6) ;
+  Int_t nSamples=TMath::Min(iBin+70,samples->GetSize()) ; //Here we set number of points to fit (70)
+  // iBin - first non-zero sample 
+  Int_t tStep=times->At(iBin+1)-times->At(iBin) ;
+  Double_t ddt=times->At(iBin)-t0-tStep ;
+  Double_t exp1=TMath::Exp(-alpha*ddt) ;
+  Double_t exp2=TMath::Exp(-beta*ddt) ;
+  Double_t dexp1=TMath::Exp(-alpha*tStep) ;
+  Double_t dexp2=TMath::Exp(-beta*tStep) ;
+  for(Int_t i = iBin; i<nSamples ; i++) {
+    Double_t dt=double(times->At(i))-t0 ;
+    Double_t fsample = double(samples->At(i)) ;
+    if(fsample>=overflow)
+      continue ;
+    Double_t diff ;
+    exp1*=dexp1 ;
+    exp2*=dexp2 ;
+    if(dt<=0.){
+      diff=fsample - ped ; 
+      fret += diff*diff ;
+      continue ;
+    }
+    Double_t dtn=TMath::Power(dt,n) ;
+    Double_t dtnE=dtn*exp1 ;
+    Double_t dt2E=dt*dt*exp2 ;
+    Double_t fit=ped+en*(dtnE + b*dt2E) ;
+    diff = fsample - fit ;
+    fret += diff*diff ;
+    if(iflag == 2){  // calculate gradient
+      Grad[0] += en*diff*(dtnE*(n/dt-alpha)+b*dt2E*(2./dt-beta))  ; //derivative over t0
+      Grad[1] -= diff*(dtnE+b*dt2E) ;
+      Grad[2] -= en*diff*dt2E ;
+    }
+  }
+  if(iflag == 2)
+    for(Int_t iparam = 0 ; iparam < 3 ; iparam++)    
+      Grad[iparam] *= 2. ; 
+}
+//-----------------------------------------------------------------------------
+Double_t AliPHOSRawFitterv1::Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * params){  //Function for fitting samples
+  //parameters:
+  //dt-time after start
+  //en-amplutude
+  //function parameters
+  
+  Double_t ped=params->At(4) ;
+  if(dt<0.)
+    return ped ; //pedestal
+  else
+    return ped+en*(TMath::Power(dt,params->At(0))*TMath::Exp(-dt*params->At(1))+b*dt*dt*TMath::Exp(-dt*params->At(3))) ;
+}
diff --git a/PHOS/AliPHOSRawFitterv1.h b/PHOS/AliPHOSRawFitterv1.h
new file mode 100644 (file)
index 0000000..f02e41d
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALIPHOSRAWFITTERV1_H
+#define ALIPHOSRAWFITTERV1_H
+/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                          */
+
+/* $Id: $ */
+
+// This class extracts the PHOS "digits" of current event
+// (amplitude,time, position,gain) from the raw stream 
+// provided by AliRawReader. See cxx source for use case.
+
+#include "AliPHOSRawFitterv0.h"
+class TArrayD;
+class TList;
+
+class AliPHOSRawFitterv1 : public AliPHOSRawFitterv0 {
+
+public:
+
+  AliPHOSRawFitterv1();
+  AliPHOSRawFitterv1(const AliPHOSRawFitterv1& rawFitter);
+  AliPHOSRawFitterv1& operator = (const AliPHOSRawFitterv1& rawFitter);
+  virtual ~AliPHOSRawFitterv1();
+
+  virtual Bool_t Eval();
+
+  static Double_t Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * fitparams) ; // Shape of correct sample
+                                                 //class member function (not object member function)
+  static void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;
+                                            // Chi^2 of the fit. Should be static to be passed to MINUIT
+  void SetLowGainParams(Int_t n, Double_t * params){fSampleParamsLow->Set(n,params) ;}  //fixed parameters of fit function
+  void SetHighGainParams(Int_t n,Double_t * params){fSampleParamsHigh->Set(n,params) ;} //fixed parameters of fit function
+
+private:
+  TArrayD *fSampleParamsLow ;   //Fixed params of sample parameterization for Low gain 
+  TArrayD *fSampleParamsHigh;   //Fixed params of sample parameterization for High gain
+  TList   *fToFit ;             //! container to transfer parameters and data to fit
+  
+  ClassDef(AliPHOSRawFitterv1,1)
+};
+
+#endif
diff --git a/PHOS/AliPHOSRawFitterv2.cxx b/PHOS/AliPHOSRawFitterv2.cxx
new file mode 100644 (file)
index 0000000..f38114c
--- /dev/null
@@ -0,0 +1,251 @@
+/**************************************************************************
+ * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+// This class plots samples and qualify their quality according to their shape
+// 
+// Typical use case:
+//     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
+//     fitter->SetSamples(sig,sigStart,sigLength);
+//     fitter->SetNBunches(nBunches);
+//     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
+//     fitter->SetCalibData(fgCalibData) ;
+//     fitter->Eval();
+//     Double_t amplitude = fitter.GetEnergy();
+//     Double_t time      = fitter.GetTime();
+//     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
+
+// Author: Dmitri Peressounko (Oct.2008)
+// Modified: Yuri Kharlov (Jul.2009)
+
+// --- ROOT system ---
+#include "TList.h"
+#include "TMath.h"
+#include "TMinuit.h"
+#include "TCanvas.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TF1.h"
+#include "TROOT.h"
+
+// --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSRawFitterv2.h"
+#include "AliPHOSPulseGenerator.h"
+
+ClassImp(AliPHOSRawFitterv2)
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv2::AliPHOSRawFitterv2():
+  AliPHOSRawFitterv0(),
+  fNtimeSamples(25),
+  fRMScut(11.)
+{
+  //Default constructor.
+  fLGpar[0]=0.971 ;
+  fLGpar[1]=0.0465;
+  fLGpar[2]=1.56  ;
+  fHGpar[0]=0.941 ; 
+  fHGpar[1]=0.0436;
+  fHGpar[2]=1.65  ;
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
+{
+  //Destructor.
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
+  AliPHOSRawFitterv0(phosFitter), 
+  fNtimeSamples(25),
+  fRMScut(11.)
+{
+  //Copy constructor.
+  fNtimeSamples=phosFitter.fNtimeSamples ;
+  for(Int_t i=0; i<3;i++){
+    fLGpar[i]=phosFitter.fLGpar[i] ;
+    fHGpar[i]=phosFitter.fHGpar[i] ;
+  }
+  fRMScut=phosFitter.fRMScut ;
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter)
+{
+  //Assignment operator.
+
+  fNtimeSamples=phosFitter.fNtimeSamples ;
+  for(Int_t i=0; i<3;i++){
+    fLGpar[i]=phosFitter.fLGpar[i] ;
+    fHGpar[i]=phosFitter.fHGpar[i] ;
+  }
+  fRMScut=phosFitter.fRMScut ;
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliPHOSRawFitterv2::Eval()
+{
+  //Extract an energy deposited in the crystal,
+  //crystal' position (module,column,row),
+  //time and gain (high or low).
+  //First collects sample, then evaluates it and if it has
+  //reasonable shape, fits it with Gamma2 function and extracts 
+  //energy and time.
+
+  if (fNBunches > 1) {
+    fQuality = 1000;
+    return kTRUE;
+  }
+
+  const Float_t kBaseLine   = 1.0;
+  const Int_t   kPreSamples = 10;
+
+  Float_t  pedMean   = 0;
+  Float_t  pedRMS    = 0;
+  Int_t    nPed      = 0;
+  UShort_t maxSample = 0;
+  Int_t    nMax      = 0;
+
+  TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ;
+  if(!cs)
+    cs = new TCanvas("CSample","CSample") ;
+
+  TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ;
+  if(!h) h = new TH1D("hSample","",200,0.,200.) ;
+
+  Double_t pedestal = 0;
+  for (Int_t i=0; i<fLength; i++) {
+    if (i<kPreSamples) {
+      nPed++;
+      pedMean += fSignal[i];
+      pedRMS  += fSignal[i]*fSignal[i] ;
+    }
+    if(fSignal[i] > maxSample) maxSample = fSignal[i];
+    if(fSignal[i] == maxSample) nMax++;
+    h->SetBinContent(i+1,fSignal[i]) ;
+  }
+  fEnergy = (Double_t)maxSample;
+  if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
+
+  if (fPedSubtract) {
+    if (nPed > 0) {
+      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
+      if(fPedestalRMS > 0.) 
+       fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
+      fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
+    }
+    else
+      return kFALSE;
+  }
+  else {
+    //take pedestals from DB
+    pedestal = (Double_t) fAmpOffset ;
+    if (fCalibData) {
+      Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
+      Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
+      pedestal += truePed - altroSettings ;
+    }
+    else{
+      AliWarning(Form("Can not read data from OCDB")) ;
+    }
+    fEnergy-=pedestal ;
+  }
+
+  if (fEnergy < kBaseLine) {
+    fEnergy = 0;
+    return kTRUE;
+  }
+  
+  // calculate time
+  fTime=0. ;
+  Double_t tRMS = 0. ;
+  Double_t tW = 0. ;
+  Int_t cnts=0 ;
+  Double_t a=0,b=0,c=0 ;
+  if(fCaloFlag == 0){ // Low gain
+    a=fLGpar[0] ; 
+    b=fLGpar[1] ; 
+    c=fLGpar[2] ; 
+  }
+  else if(fCaloFlag == 1){ // High gain
+    a=fHGpar[0] ; 
+    b=fHGpar[1] ; 
+    c=fHGpar[2] ; 
+  }
+
+
+  fQuality = 0. ;
+  
+  for(Int_t i=1; i<fLength && cnts<fNtimeSamples; i++){
+    if(fSignal[i] < pedestal)
+      continue ;
+    Double_t de = fSignal[i]   - pedestal ;
+    Double_t av = fSignal[i-1] - pedestal + de;
+    if(av<=0.) //this is fluctuation around pedestal, skip it
+      continue ;
+    Double_t ds = fSignal[i] - fSignal[i-1] ;
+    Double_t ti = ds/av ;       // calculate log. derivative
+    ti = a/(ti+b)-c*ti ;        // and compare with parameterization
+    ti = i - ti ; 
+    Double_t wi = TMath::Abs(ds) ;
+    fTime += ti*wi ;
+    tW    += wi;
+    tRMS  += ti*ti*wi ;
+    cnts++ ;
+  } 
+
+  if(tW>0.){
+    fTime/=tW ;
+    fQuality = tRMS/tW-fTime*fTime ;
+  }
+  else{
+    fTime=-999. ;
+    fQuality=999. ;
+  }
+
+  Bool_t isBad = 0 ;
+  for(Int_t i=1; i<fLength-1&&!isBad; i++){
+    if(fSignal[i] > fSignal[i-1]+5 && fSignal[i] > fSignal[i+1]+5) { //single jump
+      isBad=1 ;
+    }
+  }
+  if(pedestal < 10.)
+    isBad=1 ;
+  
+  if(fPedestalRMS > 0.1)
+    isBad=1 ;
+  
+  for(Int_t i=1; i<fLength-1&&!isBad; i++){
+    if(fSignal[i] < pedestal-1)
+      isBad=1 ;
+  }
+  
+  if(fEnergy>10. && !isBad ){
+    printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ;
+    if(fOverflow)printf(" Overflow \n") ;
+    if(isBad)printf("bad") ;
+    cs->cd() ;
+    h->Draw() ;
+    cs->Update() ;
+    getchar() ;
+  }
+
+  return kTRUE;
+}
diff --git a/PHOS/AliPHOSRawFitterv2.h b/PHOS/AliPHOSRawFitterv2.h
new file mode 100644 (file)
index 0000000..9d22362
--- /dev/null
@@ -0,0 +1,40 @@
+#ifndef ALIPHOSRAWFITTERV2_H
+#define ALIPHOSRAWFITTERV2_H
+/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                          */
+
+/* $Id: $ */
+
+// This class extracts the PHOS "digits" of current event
+// (amplitude,time, position,gain) from the raw stream 
+// provided by AliRawReader. See cxx source for use case.
+
+#include "AliPHOSRawFitterv0.h"
+class TList;
+
+class AliPHOSRawFitterv2 : public AliPHOSRawFitterv0 {
+
+public:
+
+  AliPHOSRawFitterv2();
+  AliPHOSRawFitterv2(const AliPHOSRawFitterv2& rawFitter);
+  AliPHOSRawFitterv2& operator = (const AliPHOSRawFitterv2& rawFitter);
+  virtual ~AliPHOSRawFitterv2();
+
+  virtual Bool_t Eval();
+
+  void SetNTimeSamples(Short_t n=25)     { fNtimeSamples=n ;}
+  void SetLowGainTParams (Double_t *pars){ for(Int_t i=0;i<3;i++) fLGpar[i]=pars[i] ;}
+  void SetHighGainTParams(Double_t *pars){ for(Int_t i=0;i<3;i++) fHGpar[i]=pars[i] ;}
+  void SetRMScut(Double_t cut=2.)        { fRMScut = cut ;}
+
+private:
+  Short_t  fNtimeSamples ; //Number of samples (after start) used to extract time
+  Double_t fLGpar[3] ;     //parameters for shape parameterization
+  Double_t fHGpar[3] ;     //parameters for shape parameterization
+  Double_t fRMScut ;       //cut to estmate goodness of sample
+  
+  ClassDef(AliPHOSRawFitterv2,1)
+};
+
+#endif
index c1f47db..ce326e2 100644 (file)
@@ -43,7 +43,7 @@ AliPHOSRecoParam::AliPHOSRecoParam() :
   fEMCSubtractPedestals(kTRUE),
   fEMCUnfold(kTRUE),
   fEMCEnergyCorrectionOn(kTRUE),
-  fEMCDecoderVersion(""),
+  fEMCFitterVersion(""),
   fGlobalAltroOffset(0),
   fGlobalAltroThreshold(0),
   fCPVClusteringThreshold(0.0),
@@ -69,7 +69,7 @@ AliPHOSRecoParam::AliPHOSRecoParam(const AliPHOSRecoParam& ):
   fEMCSubtractPedestals(kTRUE),
   fEMCUnfold(kTRUE),
   fEMCEnergyCorrectionOn(kTRUE),
-  fEMCDecoderVersion(""),
+  fEMCFitterVersion(""),
   fGlobalAltroOffset(0),
   fGlobalAltroThreshold(0),
   fCPVClusteringThreshold(0.0),
@@ -98,7 +98,7 @@ AliPHOSRecoParam& AliPHOSRecoParam::operator = (const AliPHOSRecoParam& recoPara
     fEMCSubtractPedestals   = recoParam.fEMCSubtractPedestals;
     fEMCUnfold              = recoParam.fEMCUnfold;
     fEMCEnergyCorrectionOn  = recoParam.fEMCEnergyCorrectionOn;
-    fEMCDecoderVersion      = recoParam.fEMCDecoderVersion;
+    fEMCFitterVersion       = recoParam.fEMCFitterVersion;
     fGlobalAltroOffset      = recoParam.fGlobalAltroOffset;
     fGlobalAltroThreshold   = recoParam.fGlobalAltroThreshold;
     fCPVClusteringThreshold = recoParam.fCPVClusteringThreshold;
@@ -126,7 +126,7 @@ void AliPHOSRecoParam::Print(Option_t * /*option*/) const
                  "\tEMCSubtractPedestals   = %d\n"
                  "\tEMCUnfold              = %d\n"
                  "\tEMCEnergyCorrectionOn  = %d\n"
-                 "\tEMCDecoderVersion      = \"%s\"\n"
+                 "\tEMCFitterVersion       = \"%s\"\n"
                  "\tGlobalAltroOffset      = %d",
                  "\tGlobalAltroThreshold   = %d",
                  fEMCClusteringThreshold,
@@ -140,7 +140,7 @@ void AliPHOSRecoParam::Print(Option_t * /*option*/) const
                  fEMCSubtractPedestals,
                  fEMCUnfold,
                  fEMCEnergyCorrectionOn,
-                 fEMCDecoderVersion.Data(),
+                 fEMCFitterVersion.Data(),
                  fGlobalAltroOffset,
                  fGlobalAltroThreshold));
 
index fb93958..bd20370 100644 (file)
@@ -29,7 +29,7 @@ public:
   Bool_t  EMCEcore2ESD()              const { return fEMCEcore2ESD;            }
   Bool_t  EMCSubtractPedestals()      const { return fEMCSubtractPedestals;    }
   Bool_t  EMCToUnfold()               const { return fEMCUnfold;               }
-  const char* EMCDecoderVersion()     const { return fEMCDecoderVersion.Data();}
+  const char* EMCFitterVersion()      const { return fEMCFitterVersion.Data(); }
   Bool_t  GetEMCEnergyCorrectionOn()  const { return fEMCEnergyCorrectionOn;   }
   Int_t   GetGlobalAltroOffset()      const { return fGlobalAltroOffset ;      }
   Int_t   GetGlobalAltroThreshold()   const { return fGlobalAltroThreshold ;   }
@@ -49,7 +49,7 @@ public:
   void SetEMCEcoreRadius(Float_t rCore)              { fEMCEcoreRadius        =rCore;   }
   void SetEMCEcore2ESD(Bool_t ecore)                 { fEMCEcore2ESD          =ecore;   }
   void SetEMCSubtractPedestals(Bool_t subtract)      { fEMCSubtractPedestals  =subtract;} 
-  void SetEMCDecoderVersion(const char* version="v1"){ fEMCDecoderVersion     =version ;}
+  void SetEMCFitterVersion(const char* version="v1") { fEMCFitterVersion     =version ; }
   void SetEMCUnfolding(Bool_t toUnfold=kFALSE)       { fEMCUnfold             =toUnfold;}
   void SetEMCEnergyCorrectionOn(Bool_t on=kTRUE)     { fEMCEnergyCorrectionOn =on;      }
   void SetGlobalAltroOffset(Int_t offset=5)          { fGlobalAltroOffset     =offset ; }
@@ -79,7 +79,7 @@ protected:
   Bool_t  fEMCSubtractPedestals;   // EMC: true if pedestal should be subtracted (in non-ZS)
   Bool_t  fEMCUnfold;              // EMC: true if overlapped clusters should be unfolded
   Bool_t  fEMCEnergyCorrectionOn;  // EMC: if true do non-linear correction of cluster energy
-  TString fEMCDecoderVersion ;     // EMC: AliPHOSRawDecoder version
+  TString fEMCFitterVersion ;      // EMC: AliPHOSRawFitter version
   Int_t   fGlobalAltroOffset ;     // Offset used in ALTRO chips in SZ runs
   Int_t   fGlobalAltroThreshold ;  // Threshold used in ALTRO chips in SZ runs
 
index c3a1f6d..fbaab3e 100644 (file)
@@ -46,9 +46,9 @@
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSEmcRecPoint.h"
 #include "AliPHOSRecParticle.h"
-#include "AliPHOSRawDecoder.h"
-#include "AliPHOSRawDecoderv1.h"
-#include "AliPHOSRawDecoderv2.h"
+#include "AliPHOSRawFitterv0.h"
+#include "AliPHOSRawFitterv1.h"
+#include "AliPHOSRawFitterv2.h"
 #include "AliPHOSRawDigiProducer.h"
 #include "AliPHOSPulseGenerator.h"
 
@@ -343,7 +343,7 @@ void  AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
   // Works on a single-event basis
   rawReader->Reset() ; 
 
-  AliPHOSRawDecoder * dc ;
+  AliPHOSRawFitterv0 * fitter ;
 
   const TObjArray* maps = AliPHOSRecoParam::GetMappings();
   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
@@ -353,55 +353,54 @@ void  AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
     mapping[i] = (AliAltroMapping*)maps->At(i);
   }
 
-  if      (strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0) 
-    dc=new AliPHOSRawDecoderv1(rawReader,mapping);
-  else if (strcmp(GetRecoParam()->EMCDecoderVersion(),"v2")==0) 
-    dc=new AliPHOSRawDecoderv2(rawReader,mapping);
+  if      (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0) 
+    fitter=new AliPHOSRawFitterv1();
+  else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0) 
+    fitter=new AliPHOSRawFitterv2();
   else
-    dc=new AliPHOSRawDecoder(rawReader,mapping);
+    fitter=new AliPHOSRawFitterv0();
 
-  dc->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
-  dc->SetAmpOffset(GetRecoParam()->GetGlobalAltroOffset());
-  dc->SetAmpThreshold(GetRecoParam()->GetGlobalAltroThreshold());
+  fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
+  fitter->SetAmpOffset     (GetRecoParam()->GetGlobalAltroOffset());
+  fitter->SetAmpThreshold  (GetRecoParam()->GetGlobalAltroThreshold());
 
   TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
   digits->SetName("DIGITS");
   Int_t bufsize = 32000;
   digitsTree->Branch("PHOS", &digits, bufsize);
 
-  //AliPHOSRawDigiProducer pr(GetRecoParam());
-  AliPHOSRawDigiProducer pr;
+  AliPHOSRawDigiProducer rdp(rawReader,mapping);
 
-  pr.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC
-  pr.SetCpvMinAmp(GetRecoParam()->GetCPVMinE());
-  pr.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut());
-  pr.MakeDigits(digits,dc);
+  rdp.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC
+  rdp.SetCpvMinAmp(GetRecoParam()->GetCPVMinE());
+  rdp.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut());
+  rdp.MakeDigits(digits,fitter);
 
-  delete dc ;
+  delete fitter ;
 
-  //!!!!for debug!!!
-/*
-  Int_t modMax=-111;
-  Int_t colMax=-111;
-  Int_t rowMax=-111;
-  Float_t eMax=-333;
-  //!!!for debug!!!
-
-  Int_t relId[4];
-  for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
-    AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
-    if(digit->GetEnergy()>eMax) {
-      fGeom->AbsToRelNumbering(digit->GetId(),relId);
-      eMax=digit->GetEnergy();
-      modMax=relId[0];
-      rowMax=relId[2];
-      colMax=relId[3];
+  if (AliLog::GetGlobalDebugLevel() == 1) {
+    Int_t modMax=-111;
+    Int_t colMax=-111;
+    Int_t rowMax=-111;
+    Float_t eMax=-333;
+    //!!!for debug!!!
+    
+    Int_t relId[4];
+    for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
+      AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
+      if(digit->GetEnergy()>eMax) {
+       fGeom->AbsToRelNumbering(digit->GetId(),relId);
+       eMax=digit->GetEnergy();
+       modMax=relId[0];
+       rowMax=relId[2];
+       colMax=relId[3];
+      }
     }
+    
+    AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
+                   modMax,colMax,rowMax,eMax));
   }
 
-  AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
-                 modMax,colMax,rowMax,eMax));
-*/
   digitsTree->Fill();
   digits->Delete();
   delete digits;
index 6aae51f..8b7f60b 100644 (file)
@@ -12,9 +12,9 @@ set(SRCS
                  AliPHOSCpvCalibData.cxx
                  AliPHOSEmcCalibData.cxx
                  AliPHOSPulseGenerator.cxx
-                 AliPHOSRawDecoder.cxx
-                 AliPHOSRawDecoderv1.cxx
-                 AliPHOSRawDecoderv2.cxx
+                 AliPHOSRawFitter.cxx
+                 AliPHOSRawFitterv1.cxx
+                 AliPHOSRawFitterv2.cxx
                  AliPHOSRawDigiProducer.cxx
                  AliPHOSEmcBadChannelsMap.cxx
                  AliPHOSSurvey.cxx
index 28f6888..b930df0 100644 (file)
@@ -18,9 +18,9 @@
 #pragma link C++ class AliPHOSCpvCalibData+;
 #pragma link C++ class AliPHOSEmcCalibData+;
 #pragma link C++ class AliPHOSPulseGenerator+;
-#pragma link C++ class AliPHOSRawDecoder+;
-#pragma link C++ class AliPHOSRawDecoderv1+;
-#pragma link C++ class AliPHOSRawDecoderv2+;
+#pragma link C++ class AliPHOSRawFitterv0+;
+#pragma link C++ class AliPHOSRawFitterv1+;
+#pragma link C++ class AliPHOSRawFitterv2+;
 #pragma link C++ class AliPHOSRawDigiProducer+;
 #pragma link C++ class AliPHOSEmcBadChannelsMap+;
 #pragma link C++ class AliPHOSSurvey+;
index af5a436..0c463c4 100644 (file)
@@ -12,9 +12,9 @@ SRCS          =  \
                  AliPHOSCpvCalibData.cxx \
                  AliPHOSEmcCalibData.cxx \
                  AliPHOSPulseGenerator.cxx \
-                 AliPHOSRawDecoder.cxx \
-                 AliPHOSRawDecoderv1.cxx \
-                 AliPHOSRawDecoderv2.cxx \
+                AliPHOSRawFitterv0.cxx \
+                AliPHOSRawFitterv1.cxx \
+                AliPHOSRawFitterv2.cxx \
                  AliPHOSRawDigiProducer.cxx \
                  AliPHOSEmcBadChannelsMap.cxx \
                  AliPHOSSurvey.cxx \