consolidate zero-length arrays (aka struct hack)
[u/mrichter/AliRoot.git] / VZERO / AliVZERO.cxx
index 1375b3e..0e69e9d 100755 (executable)
@@ -42,6 +42,8 @@
 #include "TBranch.h"
 #include "TClonesArray.h"
 #include "TStopwatch.h"
+#include "TParameter.h"
+#include "TF1.h"
 
 // --- AliRoot header files ---
 #include "AliRun.h"
 #include "AliVZEROLoader.h"
 #include "AliVZERODigitizer.h"
 #include "AliVZEROBuffer.h"
-#include "AliRunDigitizer.h"
+#include "AliDigitizationInput.h"
 #include "AliVZEROdigit.h"
+#include "AliVZEROSDigit.h"
 #include "AliDAQ.h"
 #include "AliRawReader.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
 #include "AliVZERORawStream.h"
+#include "AliVZEROCalibData.h"
+#include "AliVZERORecoParam.h"
+#include "AliVZEROReconstructor.h"
 
 ClassImp(AliVZERO)
  //__________________________________________________________________
@@ -65,13 +73,23 @@ AliVZERO::AliVZERO(): AliDetector(),
          fMaxStepQua(0.),
          fMaxStepAlu(0.),
          fMaxDestepQua(0.),
-         fMaxDestepAlu(0.)
+          fMaxDestepAlu(0.),
+          fCalibData(NULL),
+          fTimeSlewing(NULL),
+          fSignalShape(NULL),
+          fRecoParam(NULL)
 {
 /// Default Constructor
     
     AliDebug(1,Form("default (empty) ctor this=%p",this));
-    fIshunt          = 0;        
+    fIshunt          = 0;
+
+    for(Int_t i = 0 ; i < 64; ++i) {
+      fNBins[i] = 0;
+      fBinSize[i] = 0;
+    }
 }
+
 //_____________________________________________________________________________
 AliVZERO::AliVZERO(const char *name, const char *title)
        : AliDetector(name,title),
@@ -81,7 +99,11 @@ AliVZERO::AliVZERO(const char *name, const char *title)
         fMaxStepQua(0.05),
         fMaxStepAlu(0.01),
         fMaxDestepQua(-1.0),
-        fMaxDestepAlu(-1.0)
+        fMaxDestepAlu(-1.0),
+        fCalibData(NULL),
+        fTimeSlewing(NULL),
+         fSignalShape(NULL),
+        fRecoParam(NULL)
 {
   
   // Standard constructor for VZERO Detector
@@ -104,7 +126,10 @@ AliVZERO::AliVZERO(const char *name, const char *title)
 //   fMaxDestepQua =  -1.0;
 //   fMaxDestepAlu =  -1.0;
   
-  
+  for(Int_t i = 0 ; i < 64; ++i) {
+    fNBins[i] = 0;
+    fBinSize[i] = 0;
+  }
 }
 
 //_____________________________________________________________________________
@@ -123,17 +148,17 @@ AliVZERO::~AliVZERO()
         fDigits->Delete();
         delete fDigits;
         fDigits=0; }
+    if (fSignalShape) {
+      delete fSignalShape;
+      fSignalShape = NULL;
+    }
+    if (fRecoParam) {
+      delete fRecoParam;
+      fRecoParam = NULL;
+    }
 }
 
 //_____________________________________________________________________________
-void AliVZERO::BuildGeometry()
-{
-  //
-  // Builds simple ROOT TNode geometry for event display
-  //
-}
-//_____________________________________________________________________________
 void AliVZERO::CreateGeometry()
 {
   //
@@ -149,17 +174,6 @@ void AliVZERO::CreateMaterials()
 }
 
 //_____________________________________________________________________________
-Int_t AliVZERO::DistanceToPrimitive(Int_t /*px*/, Int_t /*py*/)
-{
-  //
-  // Calculates the distance from the mouse to the VZERO on the screen
-  // Dummy routine
-  //
-  
-  return 9999;
-}
-//_____________________________________________________________________________
 void AliVZERO::Init()
 {
   //
@@ -233,115 +247,176 @@ void AliVZERO::SetTreeAddress()
 }
 
 //_____________________________________________________________________________
-AliDigitizer* AliVZERO::CreateDigitizer(AliRunDigitizer* manager) const
+AliDigitizer* AliVZERO::CreateDigitizer(AliDigitizationInput* digInput) const
 {
   //
   // Creates a digitizer for VZERO
   //
-  return new AliVZERODigitizer(manager);
+  return new AliVZERODigitizer(digInput);
 }
 
 //_____________________________________________________________________________
 void AliVZERO::Hits2Digits(){
   //
-  // Converts hits to digits of the current event
+  // Converts hits to digits
   //
-  // Inputs file name
-  const char *alifile = "galice.root";   
+  // Creates the VZERO digitizer 
+  AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2Digits);
 
-  // Create the run digitizer 
-  AliRunDigitizer* manager = new AliRunDigitizer(1, 1);
-  manager->SetInputStream(0, alifile);
-  manager->SetOutputFile("H2Dfile");
+  // Creates the digits
+  dig->Digitize("");
+
+  // deletes the digitizer
+  delete dig;
+}
 
+//_____________________________________________________________________________
+void AliVZERO::Hits2SDigits(){
+  //
+  // Converts hits to summable digits
+  //
   // Creates the VZERO digitizer 
-  AliVZERODigitizer* dig = new AliVZERODigitizer(manager);
+  AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2SDigits);
 
-  // Creates the digits
-  dig->Exec("");
+  // Creates the sdigits
+  dig->Digitize("");
 
+  // deletes the digitizer
+  delete dig;
 }
+
 //_____________________________________________________________________________
 void AliVZERO::Digits2Raw()
 {
-  //
-  // Converts digits of the current event to raw data
-  //
-  AliVZERO *fVZERO = (AliVZERO*)gAlice->GetDetector("VZERO");
-  fLoader->LoadDigits();
-  TTree* digits = fLoader->TreeD();
-  if (!digits) {
-    Error("Digits2Raw", "no digits tree");
-    return;
-  }
-  TClonesArray * VZEROdigits = new TClonesArray("AliVZEROdigit",1000);
-  fVZERO->SetTreeAddress();            
-  digits->GetBranch("VZERODigit")->SetAddress(&VZEROdigits); 
+   //
+   //  Converts digits of the current event to raw data
+   //
   
-  const char *fileName    = AliDAQ::DdlFileName("VZERO",0);
-  AliVZEROBuffer* buffer  = new AliVZEROBuffer(fileName);
+   AliVZERO *fVZERO = (AliVZERO*)gAlice->GetDetector("VZERO");
+   fLoader->LoadDigits();
+   TTree* digits = fLoader->TreeD();
+   if (!digits) {
+      Error("Digits2Raw", "no digits tree");
+      return;
+   }
+   TClonesArray * VZEROdigits = new TClonesArray("AliVZEROdigit",1000);
+   fVZERO->SetTreeAddress();           
+   digits->GetBranch("VZERODigit")->SetAddress(&VZEROdigits); 
   
-  //  Verbose level
-  //  0: Silent
-  //  1: cout messages
-  //  2: txt files with digits 
-  //  BE CAREFUL, verbose level 2 MUST be used only for debugging and
-  //  it is highly suggested to use this mode only for debugging digits files
-  //  reasonably small, because otherwise the size of the txt files can reach
-  //  quickly several MB wasting time and disk space.
+   const char *fileName    = AliDAQ::DdlFileName("VZERO",0);
+   AliVZEROBuffer* buffer  = new AliVZEROBuffer(fileName);
   
-  ofstream ftxt;
-  buffer->SetVerbose(0);
-  Int_t verbose = buffer->GetVerbose();
-
-  // Get Trigger information first
-  // Read trigger inputs from trigger-detector object
-  AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
-  if( !dataLoader->IsFileOpen() ) 
-    dataLoader->OpenFile( "READ" );
-  AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get( "Trigger" );
-  UInt_t triggerInfo = 0;
-  if(trgdet) {
-    triggerInfo = trgdet->GetMask() & 0xffff;
-  }
-  else {
-    AliError(Form("There is no trigger object for %s",fLoader->GetName()));
-  }
-  buffer->WriteTriggerInfo((UInt_t)triggerInfo);
-
-  // Now write the channel information: charge+time
-  // We assume here an ordered (by PMNumber) array of
-  // digits!!
-  Int_t nEntries = Int_t(digits->GetEntries());
-  for (Int_t i = 0; i < nEntries; i++) {
+   // Get Trigger information first
+   // Read trigger inputs from trigger-detector object
+   AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
+   if( !dataLoader->IsFileOpen() ) 
+        dataLoader->OpenFile( "READ" );
+   AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get( "Trigger" );
+   UInt_t triggerInfo = 0;
+   if(trgdet) {
+      triggerInfo = trgdet->GetMask() & 0xffff;
+   }
+   else {
+      AliError(Form("There is no trigger object for %s",fLoader->GetName()));
+   }
+
+   buffer->WriteTriggerInfo((UInt_t)triggerInfo); 
+   buffer->WriteTriggerScalers(); 
+   buffer->WriteBunchNumbers(); 
   
-    fVZERO->ResetDigits();
-    digits->GetEvent(i);
-    Int_t ndig = VZEROdigits->GetEntriesFast(); 
+   Int_t aBBflagsV0A = 0;
+   Int_t aBBflagsV0C = 0;
+   Int_t aBGflagsV0A = 0;
+   Int_t aBGflagsV0C = 0;
+
+   if (digits->GetUserInfo()->FindObject("BBflagsV0A")) {
+     aBBflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0A"))->GetVal();
+   }
+   else
+     AliWarning("V0A beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
+
+   if (digits->GetUserInfo()->FindObject("BBflagsV0C")) {
+     aBBflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0C"))->GetVal();
+   }
+   else
+     AliWarning("V0C beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
+
+   if (digits->GetUserInfo()->FindObject("BGflagsV0A")) {
+     aBGflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0A"))->GetVal();
+   }
+   else
+     AliWarning("V0A beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
+
+   if (digits->GetUserInfo()->FindObject("BGflagsV0C")) {
+     aBGflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0C"))->GetVal();
+   }
+   else
+     AliWarning("V0C beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
+
+   // Now retrieve the channel information: charge smaples + time and 
+   // dump it into ADC and Time arrays
+   Int_t nEntries = Int_t(digits->GetEntries());
+   Short_t aADC[64][AliVZEROdigit::kNClocks];
+   Float_t aTime[64];
+   Float_t aWidth[64];
+   Bool_t  aIntegrator[64];
+   Bool_t  aBBflag[64];
+   Bool_t  aBGflag[64];
+  
+   for (Int_t i = 0; i < nEntries; i++) {
+     fVZERO->ResetDigits();
+     digits->GetEvent(i);
+     Int_t ndig = VZEROdigits->GetEntriesFast(); 
    
-    if(ndig == 0) continue;
-    if(verbose == 2) {ftxt.open("VZEROdigits.txt",ios::app);}
-    for(Int_t k=0; k<ndig; k++){
-        AliVZEROdigit* fVZERODigit = (AliVZEROdigit*) VZEROdigits->At(k);                      
-       Int_t ADC       = fVZERODigit->ADC();
-       Int_t PMNumber  = fVZERODigit->PMNumber();
-       Int_t Time      = fVZERODigit->Time();
-        if(verbose == 1) { cout <<"DDL: "<<fileName<< "\tdigit number: "<< k<<"\tPM number: "
-                           <<PMNumber<<"\tADC: "<< ADC << "\tTime: "<< Time << endl;} 
-       if(verbose == 2) {
-           ftxt<<"DDL: "<<fileName<< "\tdigit number: "<< k<<"\tPM number: "
-                          <<PMNumber<<"\tADC: "<< ADC << "\tTime: "<< Time << endl;          
-       }
-        buffer->WriteChannel(PMNumber, ADC, Time);
-    }
-  if(verbose==2) ftxt.close();
-  }
-
-  buffer->WriteScalers();
-  buffer->WriteMBInfo();
-
+     if(ndig == 0) continue;
+     for(Int_t k=0; k<ndig; k++){
+         AliVZEROdigit* fVZERODigit = (AliVZEROdigit*) VZEROdigits->At(k);
+         // Convert aliroot channel k into FEE channel iChannel before writing data
+        Int_t iChannel       = AliVZEROCalibData::GetBoardNumber(fVZERODigit->PMNumber()) * 8 +
+          AliVZEROCalibData::GetFEEChannelNumber(fVZERODigit->PMNumber());
+        for(Int_t iClock = 0; iClock < AliVZEROdigit::kNClocks; ++iClock) aADC[iChannel][iClock] = fVZERODigit->ChargeADC(iClock);
+        aTime[iChannel]      = fVZERODigit->Time();
+        aWidth[iChannel]     = fVZERODigit->Width();
+        aIntegrator[iChannel]= fVZERODigit->Integrator();
+        if(fVZERODigit->PMNumber() < 32) {
+          aBBflag[iChannel] = (aBBflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
+          aBGflag[iChannel] = (aBGflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
+        }
+        else {
+          aBBflag[iChannel] = (aBBflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
+          aBGflag[iChannel] = (aBGflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
+        }
+         AliDebug(1,Form("DDL: %s\tdigit number: %d\tPM number: %d\tADC: %d\tTime: %f",
+                        fileName,k,fVZERODigit->PMNumber(),aADC[iChannel][AliVZEROdigit::kNClocks/2],aTime[iChannel])); 
+     }        
+   }
+
+   // Now fill raw data
+          
+   for (Int_t  iCIU = 0; iCIU < 8; iCIU++) { 
+   // decoding of one Channel Interface Unit numbered iCIU - there are 8 channels per CIU (and 8 CIUs) :
+  
+      for(Int_t iChannel_Offset = iCIU*8; iChannel_Offset < (iCIU*8)+8; iChannel_Offset=iChannel_Offset+4) { 
+         for(Int_t iChannel = iChannel_Offset; iChannel < iChannel_Offset+4; iChannel++) {
+             buffer->WriteChannel(iChannel, aADC[iChannel], aIntegrator[iChannel]);       
+         }
+         buffer->WriteBeamFlags(&aBBflag[iChannel_Offset],&aBGflag[iChannel_Offset]); 
+         buffer->WriteMBInfo(); 
+         buffer->WriteMBFlags();   
+         buffer->WriteBeamScalers(); 
+      } 
+
+      for(Int_t iChannel = iCIU*8 + 7; iChannel >= iCIU*8; iChannel--) {
+          buffer->WriteTiming(aTime[iChannel], aWidth[iChannel]); 
+      }
+
+    // End of decoding of one CIU card
+    
+  } // end of decoding the eight CIUs
+     
   delete buffer;
-  fLoader->UnloadDigits();
+  fLoader->UnloadDigits();  
 }
 
 //_____________________________________________________________________________
@@ -354,42 +429,169 @@ Bool_t AliVZERO::Raw2SDigits(AliRawReader* rawReader){
 
   if(!fLoader) {
     AliError("no VZERO loader found");
-    return kFALSE; }
-
-  TTree* treeD  = fLoader->TreeD();
-  if(!treeD) {
-      fLoader->MakeTree("D");
-      treeD = fLoader->TreeD(); }
-        
-  AliVZEROdigit  digit;
-  AliVZEROdigit* pdigit = &digit;
-  const Int_t kBufferSize = 4000;
-   
-  treeD->Branch("VZERO", "AliVZEROdigit",  &pdigit, kBufferSize);
+    return kFALSE;
+  }
+  fLoader->LoadSDigits("UPDATE");
 
-  rawReader->Reset();
-  AliVZERORawStream* rawStream  = new AliVZERORawStream(rawReader);    
+  if (!fLoader->TreeS()) fLoader->MakeTree("S");
+  fLoader->MakeSDigitsContainer();
+  TTree* treeS  = fLoader->TreeS();
+
+  TClonesArray *sdigits = new TClonesArray("AliVZEROSDigit", 64);
+  treeS->Branch("VZEROSDigit", &sdigits); 
+
+  {
+    rawReader->Reset();
+    AliVZERORawStream rawStream(rawReader);    
      
-  if (!rawStream->Next()) return kFALSE; // No VZERO data found
-  
-  for(Int_t i=0; i<64; i++) {
-      new(pdigit) AliVZEROdigit(i, (Int_t)rawStream->GetADC(i), (Int_t)rawStream->GetTime(i)); 
-      treeD->Fill();
+    if (!rawStream.Next()) return kFALSE; // No VZERO data found
+
+    GetCalibData();
+
+    Int_t nSDigits = 0;
+    Float_t *charges = NULL;
+    Int_t nbins = 0;
+    for(Int_t iChannel=0; iChannel < 64; ++iChannel) {
+      Int_t offlineCh = rawStream.GetOfflineChannel(iChannel);
+      Short_t chargeADC[AliVZEROdigit::kNClocks];
+      for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
+       chargeADC[iClock] = rawStream.GetPedestal(iChannel,iClock);
+      }
+      // Integrator flag
+      Bool_t integrator = rawStream.GetIntegratorFlag(iChannel,AliVZEROdigit::kNClocks/2);
+      // HPTDC data (leading time and width)
+      Int_t board = AliVZEROCalibData::GetBoardNumber(offlineCh);
+      Float_t time = rawStream.GetTime(iChannel)*fCalibData->GetTimeResolution(board);
+      //      Float_t width = rawStream.GetWidth(iChannel)*fCalibData->GetWidthResolution(board);
+      Float_t adc = 0;
+
+      // Pedestal retrieval and suppression
+      Float_t maxadc = 0;
+      Int_t imax = -1;
+      Float_t adcPedSub[AliVZEROdigit::kNClocks];
+      Float_t integral = fSignalShape->Integral(0,200);
+      for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
+       Bool_t iIntegrator = (iClock%2 == 0) ? integrator : !integrator;
+       Int_t k = offlineCh + 64*iIntegrator;
+       adcPedSub[iClock] = (Float_t)chargeADC[iClock] - fCalibData->GetPedestal(k);
+       if(adcPedSub[iClock] <= fRecoParam->GetNSigmaPed()*fCalibData->GetSigma(k)) {
+         adcPedSub[iClock] = 0;
+         continue;
+       }
+       if(iClock < fRecoParam->GetStartClock() || iClock > fRecoParam->GetEndClock()) continue;
+       if(adcPedSub[iClock] > maxadc) {
+         maxadc = adcPedSub[iClock];
+         imax   = iClock;
+       }
+      }
+      if (imax != -1) {
+       Int_t start = imax - fRecoParam->GetNPreClocks();
+       if (start < 0) start = 0;
+       Int_t end = imax + fRecoParam->GetNPostClocks();
+       if (end > 20) end = 20;
+       for(Int_t iClock = start; iClock <= end; iClock++) {
+         adc += adcPedSub[iClock];
+       }
+      }
+      Float_t correctedTime = CorrectLeadingTime(offlineCh,time,adc);
+
+      if (!charges) {
+       nbins = fNBins[offlineCh];
+       charges = new Float_t[nbins];
+      }
+      else if (nbins != fNBins[offlineCh]) {
+       delete [] charges;
+       nbins = fNBins[offlineCh];
+       charges = new Float_t[nbins];
+      }
+      memset(charges,0,nbins*sizeof(Float_t));
+
+      // Now lets produce SDigit
+      if ((correctedTime > (AliVZEROReconstructor::kInvalidTime + 1e-6)) &&
+         (adc > 1e-6)) {
+       for(Int_t iBin = 0; iBin < nbins; ++iBin) {
+         Float_t t = fBinSize[offlineCh]*Float_t(iBin);
+         if ((t < correctedTime) ||
+             (t > (correctedTime+200.))) continue;
+         charges[iBin] = kChargePerADC*adc*(fSignalShape->Eval(t-correctedTime)*fBinSize[offlineCh]/integral);
+       }
+      }
+
+      TClonesArray &sdigitsref = *sdigits;  
+      new (sdigitsref[nSDigits++]) AliVZEROSDigit(offlineCh,fNBins[offlineCh],charges);
+    }
+    if (charges) delete [] charges;
   }
-// Checks if everything is OK by printing results 
-
-//   for(int i=0;i<64;i++) {
-//     printf("Channel %d : %d %d \n",i,rawStream->GetADC(i),rawStream->GetTime(i)); }
-//   treeD->Print(); printf(" \n"); 
-       
-  fLoader->WriteDigits("OVERWRITE");
-  fLoader->UnloadDigits();     
+  
+  treeS->Fill();
+  fLoader->WriteSDigits("OVERWRITE");
+  fLoader->UnloadSDigits();    
        
-  delete rawStream;
-
   timer.Stop();
   timer.Print();
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AliVZERO::GetCalibData()
+{
+  // Gets calibration object for VZERO set
+  // Do nothing in case it is already loaded
+  if (fCalibData) return;
+
+  AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/Data");
+  if (entry) fCalibData = (AliVZEROCalibData*) entry->GetObject();
+  if (!fCalibData)  AliFatal("No calibration data from calibration database !");
+
+  AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeSlewing");
+  if (!entry2) AliFatal("VZERO time slewing function is not found in OCDB !");
+  fTimeSlewing = (TF1*)entry2->GetObject();
+
+  for(Int_t i = 0 ; i < 64; ++i) {
+    Int_t board = AliVZEROCalibData::GetBoardNumber(i);
+    fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
+                            (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/
+                           fCalibData->GetTimeResolution(board));
+    fBinSize[i] = fCalibData->GetTimeResolution(board);
+  }
+
+  fSignalShape = new TF1("VZEROSDigitSignalShape",this,&AliVZERO::SignalShape,0,200,6,"AliVZERO","SignalShape");
+  fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,
+                             2.9,6.40982,3.69339e-01);
+
+  fRecoParam = new AliVZERORecoParam;
+
+  return;
 }
 
+Float_t AliVZERO::CorrectLeadingTime(Int_t i, Float_t time, Float_t adc) const
+{
+  // Correct the leading time
+  // for slewing effect and
+  // misalignment of the channels
+  if (time < 1e-6) return -1024;
+
+  // In case of pathological signals
+  if (adc < 1e-6) return time;
+
+  // Slewing correction
+  Float_t thr = fCalibData->GetCalibDiscriThr(i,kTRUE);
+  time -= fTimeSlewing->Eval(adc/thr);
 
+  return time;
+}
+
+double AliVZERO::SignalShape(double *x, double *par)
+{
+  // this function simulates the signal
+  // shape used in Raw->SDigits method
+
+  Double_t xx = x[0];
+  if (xx <= par[0]) return 0;
+  Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
+  if (xx <= par[3]) return a;
+  Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
+  Double_t f = a*b/(a+b);
+  AliDebug(100,Form("x=%f func=%f",xx,f));
+  return f;
+}