VZERO/AliVZERO.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 //                                                                       //
  20 //                          V-Zero   Detector                            //
  21 //  This class contains the base procedures for the VZERO  detector      //
  22 //  Default geometry of November 2003 :   V0R box is 4.4 cm thick        //
  23 //                                  scintillators are 2 cm thick         //
  24 //  All comments should be sent to Brigitte CHEYNIS :                    //
  25 //                                 b.cheynis@ipnl.in2p3.fr               //
  26 //                                                                       //
  27 //                                                                       //
  28 ///////////////////////////////////////////////////////////////////////////
  29
  30
  31 // --- Standard libraries ---
  32 #include <Riostream.h>
  33 #include <stdlib.h>
  34
  35 // --- ROOT libraries ---
  36 #include <TNamed.h>
  37 #include "TROOT.h"
  38 #include "TFile.h"
  39 #include "TNetFile.h"
  40 #include "TRandom.h"
  41 #include "TTree.h"
  42 #include "TBranch.h"
  43 #include "TClonesArray.h"
  44 #include "TStopwatch.h"
  45 #include "TParameter.h"
  46 #include "TF1.h"
  47
  48 // --- AliRoot header files ---
  49 #include "AliRun.h"
  50 #include "AliMC.h"
  51 #include "AliVZERO.h"
  52 #include "AliVZEROLoader.h"
  53 #include "AliVZERODigitizer.h"
  54 #include "AliVZEROBuffer.h"
  55 #include "AliRunDigitizer.h"
  56 #include "AliVZEROdigit.h"
  57 #include "AliVZEROSDigit.h"
  58 #include "AliDAQ.h"
  59 #include "AliRawReader.h"
  60 #include "AliCDBManager.h"
  61 #include "AliCDBEntry.h"
  62 #include "AliVZERORawStream.h"
  63 #include "AliVZEROCalibData.h"
  64 #include "AliVZERORecoParam.h"
  65 #include "AliVZEROReconstructor.h"
  66
  67 ClassImp(AliVZERO)
  68  //__________________________________________________________________
  69 AliVZERO::AliVZERO(): AliDetector(),
  70           fIdSens1(0),
  71           fThickness(0.),
  72           fThickness1(0.),
  73           fMaxStepQua(0.),
  74           fMaxStepAlu(0.),
  75           fMaxDestepQua(0.),
  76           fMaxDestepAlu(0.),
  77           fCalibData(NULL),
  78           fTimeSlewing(NULL),
  79           fSignalShape(NULL),
  80           fRecoParam(NULL)
  81 {
  82 /// Default Constructor
  83
  84     AliDebug(1,Form("default (empty) ctor this=%p",this));
  85     fIshunt          = 0;
  86
  87     for(Int_t i = 0 ; i < 64; ++i) {
  88       fNBins[i] = 0;
  89       fBinSize[i] = 0;
  90     }
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliVZERO::AliVZERO(const char *name, const char *title)
  95        : AliDetector(name,title),
  96          fIdSens1(0),
  97          fThickness(4.4),
  98          fThickness1(2.0),
  99          fMaxStepQua(0.05),
 100          fMaxStepAlu(0.01),
 101          fMaxDestepQua(-1.0),
 102          fMaxDestepAlu(-1.0),
 103          fCalibData(NULL),
 104          fTimeSlewing(NULL),
 105          fSignalShape(NULL),
 106          fRecoParam(NULL)
 107 {
 108
 109   // Standard constructor for VZERO Detector
 110
 111   AliDebug(1,Form("ctor this=%p",this));
 112
 113   //  fIshunt       =  1;  // All hits are associated with primary particles
 114
 115   fHits         =  new TClonesArray("AliVZEROhit", 400);
 116   fDigits       =  new TClonesArray("AliVZEROdigit",400);
 117
 118   gAlice->GetMCApp()->AddHitList(fHits);
 119
 120 //   fThickness    =  4.4;   // total thickness of the V0R box in cm
 121 //   fThickness1   =  2.0;   // thickness of scintillating cells in cm
 122 //
 123 //   fMaxStepQua   =  0.05;
 124 //   fMaxStepAlu   =  0.01;
 125 //
 126 //   fMaxDestepQua =  -1.0;
 127 //   fMaxDestepAlu =  -1.0;
 128
 129   for(Int_t i = 0 ; i < 64; ++i) {
 130     fNBins[i] = 0;
 131     fBinSize[i] = 0;
 132   }
 133 }
 134
 135 //_____________________________________________________________________________
 136 AliVZERO::~AliVZERO()
 137 {
 138   //
 139   // Default destructor for VZERO Detector
 140   //
 141
 142     if (fHits) {
 143         fHits->Delete();
 144         delete fHits;
 145         fHits=0; }
 146
 147     if (fDigits) {
 148         fDigits->Delete();
 149         delete fDigits;
 150         fDigits=0; }
 151     if (fSignalShape) {
 152       delete fSignalShape;
 153       fSignalShape = NULL;
 154     }
 155     if (fRecoParam) {
 156       delete fRecoParam;
 157       fRecoParam = NULL;
 158     }
 159 }
 160
 161 //_____________________________________________________________________________
 162 void AliVZERO::CreateGeometry()
 163 {
 164   //
 165   // Builds simple Geant3 geometry
 166   //
 167 }
 168 //_____________________________________________________________________________
 169 void AliVZERO::CreateMaterials()
 170 {
 171   //
 172   // Creates materials used for Geant3 geometry
 173   //
 174 }
 175
 176 //_____________________________________________________________________________
 177 void AliVZERO::Init()
 178 {
 179   //
 180   // Initialises the VZERO  class after it has been built
 181   //
 182 }
 183
 184
 185 //_____________________________________________________________________________
 186 void AliVZERO::SetMaxStepQua(Float_t p1)
 187 {
 188   //
 189   // Possible parametrisation of steps in active materials
 190   //
 191      fMaxStepQua = p1;
 192 }
 193
 194 //_____________________________________________________________________________
 195 void AliVZERO::SetMaxStepAlu(Float_t p1)
 196 {
 197   //
 198   // Possible parametrisation of steps in Aluminum foils (not used in
 199   // version v2)
 200   //
 201     fMaxStepAlu = p1;
 202 }
 203
 204 //_____________________________________________________________________________
 205 void AliVZERO::SetMaxDestepQua(Float_t p1)
 206 {
 207   //
 208   // Possible parametrisation of steps in active materials (quartz)
 209   //
 210     fMaxDestepQua = p1;
 211 }
 212
 213 //_____________________________________________________________________________
 214 void AliVZERO::SetMaxDestepAlu(Float_t p1)
 215 {
 216   //
 217   // Possible parametrisation of steps in Aluminum (not used in
 218   // version v2)
 219   //
 220     fMaxDestepAlu = p1;
 221 }
 222
 223 //_____________________________________________________________________________
 224 AliLoader* AliVZERO::MakeLoader(const char* topfoldername)
 225 {
 226   //
 227   // Builds VZEROgetter (AliLoader type)
 228   // if detector wants to use customized getter, it must overload this method
 229   //
 230 //  Info("MakeLoader","Creating AliVZEROLoader. Top folder is %s.",topfoldername);
 231
 232   AliDebug(1,Form("Creating AliVZEROLoader, Top folder is %s ",topfoldername));
 233   fLoader = new AliVZEROLoader(GetName(),topfoldername);
 234   return fLoader;
 235 }
 236
 237 //_____________________________________________________________________________
 238 void AliVZERO::SetTreeAddress()
 239 {
 240   //
 241   // Sets tree address for hits.
 242   //
 243   if (fLoader->TreeH() && (fHits == 0x0))
 244     fHits = new  TClonesArray("AliVZEROhit", 400);
 245
 246   AliDetector::SetTreeAddress();
 247 }
 248
 249 //_____________________________________________________________________________
 250 AliDigitizer* AliVZERO::CreateDigitizer(AliRunDigitizer* manager) const
 251 {
 252   //
 253   // Creates a digitizer for VZERO
 254   //
 255   return new AliVZERODigitizer(manager);
 256 }
 257
 258 //_____________________________________________________________________________
 259 void AliVZERO::Hits2Digits(){
 260   //
 261   // Converts hits to digits
 262   //
 263   // Creates the VZERO digitizer
 264   AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2Digits);
 265
 266   // Creates the digits
 267   dig->Exec("");
 268
 269 }
 270
 271 //_____________________________________________________________________________
 272 void AliVZERO::Hits2SDigits(){
 273   //
 274   // Converts hits to summable digits
 275   //
 276   // Creates the VZERO digitizer
 277   AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2SDigits);
 278
 279   // Creates the sdigits
 280   dig->Exec("");
 281
 282 }
 283
 284 //_____________________________________________________________________________
 285 void AliVZERO::Digits2Raw()
 286 {
 287    //
 288    //  Converts digits of the current event to raw data
 289    //
 290
 291    AliVZERO *fVZERO = (AliVZERO*)gAlice->GetDetector("VZERO");
 292    fLoader->LoadDigits();
 293    TTree* digits = fLoader->TreeD();
 294    if (!digits) {
 295       Error("Digits2Raw", "no digits tree");
 296       return;
 297    }
 298    TClonesArray * VZEROdigits = new TClonesArray("AliVZEROdigit",1000);
 299    fVZERO->SetTreeAddress();
 300    digits->GetBranch("VZERODigit")->SetAddress(&VZEROdigits);
 301
 302    const char *fileName    = AliDAQ::DdlFileName("VZERO",0);
 303    AliVZEROBuffer* buffer  = new AliVZEROBuffer(fileName);
 304
 305    // Get Trigger information first
 306    // Read trigger inputs from trigger-detector object
 307    AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
 308    if( !dataLoader->IsFileOpen() )
 309         dataLoader->OpenFile( "READ" );
 310    AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get( "Trigger" );
 311    UInt_t triggerInfo = 0;
 312    if(trgdet) {
 313       triggerInfo = trgdet->GetMask() & 0xffff;
 314    }
 315    else {
 316       AliError(Form("There is no trigger object for %s",fLoader->GetName()));
 317    }
 318
 319    buffer->WriteTriggerInfo((UInt_t)triggerInfo);
 320    buffer->WriteTriggerScalers();
 321    buffer->WriteBunchNumbers();
 322
 323    Int_t aBBflagsV0A = 0;
 324    Int_t aBBflagsV0C = 0;
 325    Int_t aBGflagsV0A = 0;
 326    Int_t aBGflagsV0C = 0;
 327
 328    if (digits->GetUserInfo()->FindObject("BBflagsV0A")) {
 329      aBBflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0A"))->GetVal();
 330    }
 331    else
 332      AliWarning("V0A beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
 333
 334    if (digits->GetUserInfo()->FindObject("BBflagsV0C")) {
 335      aBBflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0C"))->GetVal();
 336    }
 337    else
 338      AliWarning("V0C beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
 339
 340    if (digits->GetUserInfo()->FindObject("BGflagsV0A")) {
 341      aBGflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0A"))->GetVal();
 342    }
 343    else
 344      AliWarning("V0A beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
 345
 346    if (digits->GetUserInfo()->FindObject("BGflagsV0C")) {
 347      aBGflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0C"))->GetVal();
 348    }
 349    else
 350      AliWarning("V0C beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
 351
 352    // Now retrieve the channel information: charge smaples + time and
 353    // dump it into ADC and Time arrays
 354    Int_t nEntries = Int_t(digits->GetEntries());
 355    Short_t aADC[64][AliVZEROdigit::kNClocks];
 356    Float_t aTime[64];
 357    Float_t aWidth[64];
 358    Bool_t  aIntegrator[64];
 359    Bool_t  aBBflag[64];
 360    Bool_t  aBGflag[64];
 361
 362    for (Int_t i = 0; i < nEntries; i++) {
 363      fVZERO->ResetDigits();
 364      digits->GetEvent(i);
 365      Int_t ndig = VZEROdigits->GetEntriesFast();
 366
 367      if(ndig == 0) continue;
 368      for(Int_t k=0; k<ndig; k++){
 369          AliVZEROdigit* fVZERODigit = (AliVZEROdigit*) VZEROdigits->At(k);
 370          // Convert aliroot channel k into FEE channel iChannel before writing data
 371          Int_t iChannel       = AliVZEROCalibData::GetBoardNumber(fVZERODigit->PMNumber()) * 8 +
 372            AliVZEROCalibData::GetFEEChannelNumber(fVZERODigit->PMNumber());
 373          for(Int_t iClock = 0; iClock < AliVZEROdigit::kNClocks; ++iClock) aADC[iChannel][iClock] = fVZERODigit->ChargeADC(iClock);
 374          aTime[iChannel]      = fVZERODigit->Time();
 375          aWidth[iChannel]     = fVZERODigit->Width();
 376          aIntegrator[iChannel]= fVZERODigit->Integrator();
 377          if(fVZERODigit->PMNumber() < 32) {
 378            aBBflag[iChannel] = (aBBflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
 379            aBGflag[iChannel] = (aBGflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
 380          }
 381          else {
 382            aBBflag[iChannel] = (aBBflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
 383            aBGflag[iChannel] = (aBGflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
 384          }
 385          AliDebug(1,Form("DDL: %s\tdigit number: %d\tPM number: %d\tADC: %d\tTime: %f",
 386                          fileName,k,fVZERODigit->PMNumber(),aADC[k][AliVZEROdigit::kNClocks/2],aTime[k]));
 387      }
 388    }
 389
 390    // Now fill raw data
 391
 392    for (Int_t  iCIU = 0; iCIU < 8; iCIU++) {
 393
 394    // decoding of one Channel Interface Unit numbered iCIU - there are 8 channels per CIU (and 8 CIUs) :
 395
 396       for(Int_t iChannel_Offset = iCIU*8; iChannel_Offset < (iCIU*8)+8; iChannel_Offset=iChannel_Offset+4) {
 397          for(Int_t iChannel = iChannel_Offset; iChannel < iChannel_Offset+4; iChannel++) {
 398              buffer->WriteChannel(iChannel, aADC[iChannel], aIntegrator[iChannel]);
 399          }
 400          buffer->WriteBeamFlags(&aBBflag[iChannel_Offset],&aBGflag[iChannel_Offset]);
 401          buffer->WriteMBInfo();
 402          buffer->WriteMBFlags();
 403          buffer->WriteBeamScalers();
 404       }
 405
 406       for(Int_t iChannel = iCIU*8 + 7; iChannel >= iCIU*8; iChannel--) {
 407           buffer->WriteTiming(aTime[iChannel], aWidth[iChannel]);
 408       }
 409
 410     // End of decoding of one CIU card
 411
 412   } // end of decoding the eight CIUs
 413
 414   delete buffer;
 415   fLoader->UnloadDigits();
 416 }
 417
 418 //_____________________________________________________________________________
 419 Bool_t AliVZERO::Raw2SDigits(AliRawReader* rawReader){
 420   // Converts the VZERO raw data into digits
 421   // The method is used for merging simulated and
 422   // real data events
 423   TStopwatch timer;
 424   timer.Start();
 425
 426   if(!fLoader) {
 427     AliError("no VZERO loader found");
 428     return kFALSE;
 429   }
 430   fLoader->LoadSDigits("UPDATE");
 431
 432   if (!fLoader->TreeS()) fLoader->MakeTree("S");
 433   fLoader->MakeSDigitsContainer();
 434   TTree* treeS  = fLoader->TreeS();
 435
 436   TClonesArray *sdigits = new TClonesArray("AliVZEROSDigit", 64);
 437   treeS->Branch("VZEROSDigit", &sdigits);
 438
 439   {
 440     rawReader->Reset();
 441     AliVZERORawStream rawStream(rawReader);
 442
 443     if (!rawStream.Next()) return kFALSE; // No VZERO data found
 444
 445     GetCalibData();
 446
 447     Int_t nSDigits = 0;
 448     Float_t *charges = NULL;
 449     Int_t nbins = 0;
 450     for(Int_t iChannel=0; iChannel < 64; ++iChannel) {
 451       Int_t offlineCh = rawStream.GetOfflineChannel(iChannel);
 452       Short_t chargeADC[AliVZEROdigit::kNClocks];
 453       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
 454         chargeADC[iClock] = rawStream.GetPedestal(iChannel,iClock);
 455       }
 456       // Integrator flag
 457       Bool_t integrator = rawStream.GetIntegratorFlag(iChannel,AliVZEROdigit::kNClocks/2);
 458       // HPTDC data (leading time and width)
 459       Int_t board = AliVZEROCalibData::GetBoardNumber(offlineCh);
 460       Float_t time = rawStream.GetTime(iChannel)*fCalibData->GetTimeResolution(board);
 461       //      Float_t width = rawStream.GetWidth(iChannel)*fCalibData->GetWidthResolution(board);
 462       Float_t adc = 0;
 463
 464       // Pedestal retrieval and suppression
 465       Float_t maxadc = 0;
 466       Int_t imax = -1;
 467       Float_t adcPedSub[AliVZEROdigit::kNClocks];
 468       Float_t integral = fSignalShape->Integral(0,200);
 469       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
 470         Bool_t iIntegrator = (iClock%2 == 0) ? integrator : !integrator;
 471         Int_t k = offlineCh + 64*iIntegrator;
 472         adcPedSub[iClock] = (Float_t)chargeADC[iClock] - fCalibData->GetPedestal(k);
 473         if(adcPedSub[iClock] <= fRecoParam->GetNSigmaPed()*fCalibData->GetSigma(k)) {
 474           adcPedSub[iClock] = 0;
 475           continue;
 476         }
 477         if(iClock < fRecoParam->GetStartClock() || iClock > fRecoParam->GetEndClock()) continue;
 478         if(adcPedSub[iClock] > maxadc) {
 479           maxadc = adcPedSub[iClock];
 480           imax   = iClock;
 481         }
 482       }
 483       if (imax != -1) {
 484         Int_t start = imax - fRecoParam->GetNPreClocks();
 485         if (start < 0) start = 0;
 486         Int_t end = imax + fRecoParam->GetNPostClocks();
 487         if (end > 20) end = 20;
 488         for(Int_t iClock = start; iClock <= end; iClock++) {
 489           adc += adcPedSub[iClock];
 490         }
 491       }
 492       Float_t correctedTime = CorrectLeadingTime(offlineCh,time,adc);
 493
 494       if (!charges) {
 495         nbins = fNBins[offlineCh];
 496         charges = new Float_t[nbins];
 497       }
 498       else if (nbins != fNBins[offlineCh]) {
 499         delete [] charges;
 500         nbins = fNBins[offlineCh];
 501         charges = new Float_t[nbins];
 502       }
 503       memset(charges,0,nbins*sizeof(Float_t));
 504
 505       // Now lets produce SDigit
 506       if ((correctedTime > (AliVZEROReconstructor::kInvalidTime + 1e-6)) &&
 507           (adc > 1e-6)) {
 508         for(Int_t iBin = 0; iBin < nbins; ++iBin) {
 509           Float_t t = fBinSize[offlineCh]*Float_t(iBin);
 510           if ((t < correctedTime) ||
 511               (t > (correctedTime+200.))) continue;
 512           charges[iBin] = kChargePerADC*adc*(fSignalShape->Eval(t-correctedTime)*fBinSize[offlineCh]/integral);
 513         }
 514       }
 515
 516       TClonesArray &sdigitsref = *sdigits;
 517       new (sdigitsref[nSDigits++]) AliVZEROSDigit(offlineCh,fNBins[offlineCh],charges);
 518     }
 519     if (charges) delete [] charges;
 520   }
 521
 522   treeS->Fill();
 523   fLoader->WriteSDigits("OVERWRITE");
 524   fLoader->UnloadSDigits();
 525
 526   timer.Stop();
 527   timer.Print();
 528   return kTRUE;
 529 }
 530
 531 //_____________________________________________________________________________
 532 void AliVZERO::GetCalibData()
 533 {
 534   // Gets calibration object for VZERO set
 535   // Do nothing in case it is already loaded
 536   if (fCalibData) return;
 537
 538   AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/Data");
 539   if (entry) fCalibData = (AliVZEROCalibData*) entry->GetObject();
 540   if (!fCalibData)  AliFatal("No calibration data from calibration database !");
 541
 542   AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeSlewing");
 543   if (!entry2) AliFatal("VZERO time slewing function is not found in OCDB !");
 544   fTimeSlewing = (TF1*)entry2->GetObject();
 545
 546   for(Int_t i = 0 ; i < 64; ++i) {
 547     Int_t board = AliVZEROCalibData::GetBoardNumber(i);
 548     fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
 549                              (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/
 550                             fCalibData->GetTimeResolution(board));
 551     fBinSize[i] = fCalibData->GetTimeResolution(board);
 552   }
 553
 554   fSignalShape = new TF1("VZEROSDigitSignalShape",this,&AliVZERO::SignalShape,0,200,6,"AliVZERO","SignalShape");
 555   fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,
 556                               2.9,6.40982,3.69339e-01);
 557
 558   fRecoParam = new AliVZERORecoParam;
 559
 560   return;
 561 }
 562
 563 Float_t AliVZERO::CorrectLeadingTime(Int_t i, Float_t time, Float_t adc) const
 564 {
 565   // Correct the leading time
 566   // for slewing effect and
 567   // misalignment of the channels
 568   if (time < 1e-6) return -1024;
 569
 570   // In case of pathological signals
 571   if (adc < 1e-6) return time;
 572
 573   // Slewing correction
 574   Float_t thr = fCalibData->GetCalibDiscriThr(i,kTRUE);
 575   time -= fTimeSlewing->Eval(adc/thr);
 576
 577   return time;
 578 }
 579
 580 double AliVZERO::SignalShape(double *x, double *par)
 581 {
 582   // this function simulates the signal
 583   // shape used in Raw->SDigits method
 584
 585   Double_t xx = x[0];
 586   if (xx <= par[0]) return 0;
 587   Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
 588   if (xx <= par[3]) return a;
 589   Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
 590   Double_t f = a*b/(a+b);
 591   AliDebug(100,Form("x=%f func=%f",xx,f));
 592   return f;
 593 }