PHOS/AliPHOSClusterizerv1.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 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.97  2006/08/29 11:41:19  kharlov
  22  * Missing implementation of ctors and = operator are added
  23  *
  24  * Revision 1.96  2006/08/25 16:56:30  kharlov
  25  * Compliance with Effective C++
  26  *
  27  * Revision 1.95  2006/08/11 12:36:26  cvetan
  28  * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
  29  *
  30  * Revision 1.94  2006/08/07 12:27:49  hristov
  31  * Removing obsolete code which affected the event numbering scheme
  32  *
  33  * Revision 1.93  2006/08/01 12:20:17  cvetan
  34  * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
  35  *
  36  * Revision 1.92  2006/04/29 20:26:46  hristov
  37  * Separate EMC and CPV calibration (Yu.Kharlov)
  38  *
  39  * Revision 1.91  2006/04/22 10:30:17  hristov
  40  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
  41  *
  42  * Revision 1.90  2006/04/11 15:22:59  hristov
  43  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
  44  *
  45  * Revision 1.89  2006/03/13 14:05:42  kharlov
  46  * Calibration objects for EMC and CPV
  47  *
  48  * Revision 1.88  2006/01/11 08:54:52  hristov
  49  * Additional protection in case no calibration entry was found
  50  *
  51  * Revision 1.87  2005/11/22 08:46:43  kharlov
  52  * Updated to new CDB (Boris Polichtchouk)
  53  *
  54  * Revision 1.86  2005/11/14 21:52:43  hristov
  55  * Coding conventions
  56  *
  57  * Revision 1.85  2005/09/27 16:08:08  hristov
  58  * New version of CDB storage framework (A.Colla)
  59  *
  60  * Revision 1.84  2005/09/21 10:02:47  kharlov
  61  * Reading calibration from CDB (Boris Polichtchouk)
  62  *
  63  * Revision 1.82  2005/09/02 15:43:13  kharlov
  64  * Add comments in GetCalibrationParameters and Calibrate
  65  *
  66  * Revision 1.81  2005/09/02 14:32:07  kharlov
  67  * Calibration of raw data
  68  *
  69  * Revision 1.80  2005/08/24 15:31:36  kharlov
  70  * Setting raw digits flag
  71  *
  72  * Revision 1.79  2005/07/25 15:53:53  kharlov
  73  * Read raw data
  74  *
  75  * Revision 1.78  2005/05/28 14:19:04  schutz
  76  * Compilation warnings fixed by T.P.
  77  *
  78  */
  79
  80 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  81 //////////////////////////////////////////////////////////////////////////////
  82 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  83 //  unfolds the clusters having several local maxima.
  84 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
  85 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
  86 //  parameters including input digits branch title, thresholds etc.)
  87 //  This TTask is normally called from Reconstructioner, but can as well be used in
  88 //  standalone mode.
  89 // Use Case:
  90 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
  91 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  92 //               // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
  93 //               // and saves recpoints in branch named "recpointsname" (default = "digitsname")
  94 //  root [1] cl->ExecuteTask()
  95 //               //finds RecPoints in all events stored in galice.root
  96 //  root [2] cl->SetDigitsBranch("digits2")
  97 //               //sets another title for Digitis (input) branch
  98 //  root [3] cl->SetRecPointsBranch("recp2")
  99 //               //sets another title four output branches
 100 //  root [4] cl->SetEmcLocalMaxCut(0.03)
 101 //               //set clusterization parameters
 102 //  root [5] cl->ExecuteTask("deb all time")
 103 //               //once more finds RecPoints options are
 104 //               // deb - print number of found rec points
 105 //               // deb all - print number of found RecPoints and some their characteristics
 106 //               // time - print benchmarking results
 107
 108 // --- ROOT system ---
 109
 110 #include "TMath.h"
 111 #include "TMinuit.h"
 112 #include "TTree.h"
 113 #include "TBenchmark.h"
 114
 115 // --- Standard library ---
 116
 117 // --- AliRoot header files ---
 118 #include "AliLog.h"
 119 #include "AliPHOSGetter.h"
 120 #include "AliPHOSGeometry.h"
 121 #include "AliPHOSClusterizerv1.h"
 122 #include "AliPHOSEmcRecPoint.h"
 123 #include "AliPHOSCpvRecPoint.h"
 124 #include "AliPHOSDigit.h"
 125 #include "AliPHOSDigitizer.h"
 126 #include "AliPHOSCalibrationDB.h"
 127 #include "AliCDBManager.h"
 128 #include "AliCDBStorage.h"
 129 #include "AliCDBEntry.h"
 130
 131 ClassImp(AliPHOSClusterizerv1)
 132
 133 //____________________________________________________________________________
 134 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
 135   AliPHOSClusterizer(),
 136   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 137   fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
 138   fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
 139   fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
 140   fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
 141   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 142   fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
 143   fIsOldRCUFormat(0)
 144 {
 145   // default ctor (to be used mainly by Streamer)
 146
 147   InitParameters() ;
 148   fDefaultInit = kTRUE ;
 149 }
 150
 151 //____________________________________________________________________________
 152 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
 153   AliPHOSClusterizer(alirunFileName, eventFolderName),
 154   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 155   fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
 156   fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
 157   fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
 158   fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
 159   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 160   fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
 161   fIsOldRCUFormat(0)
 162 {
 163   // ctor with the indication of the file where header Tree and digits Tree are stored
 164
 165   InitParameters() ;
 166   Init() ;
 167   fDefaultInit = kFALSE ;
 168 }
 169
 170 //____________________________________________________________________________
 171 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
 172   AliPHOSClusterizer(obj),
 173   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 174   fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
 175   fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
 176   fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
 177   fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
 178   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 179   fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
 180   fIsOldRCUFormat(0)
 181 {
 182   // Copy constructor
 183 }
 184 //____________________________________________________________________________
 185   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 186 {
 187   // dtor
 188
 189 }
 190 //____________________________________________________________________________
 191 const TString AliPHOSClusterizerv1::BranchName() const
 192 {
 193   return GetName();
 194 }
 195
 196 //____________________________________________________________________________
 197 Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
 198 {
 199   // Convert EMC measured amplitude into real energy.
 200   // Calibration parameters are taken from calibration data base for raw data,
 201   // or from digitizer parameters for simulated data.
 202
 203   if(fCalibData){
 204     Int_t relId[4];
 205     AliPHOSGetter *gime = AliPHOSGetter::Instance();
 206     gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
 207     Int_t   module = relId[0];
 208     Int_t   column = relId[3];
 209     Int_t   row    = relId[2];
 210     if(absId <= fEmcCrystals) { // this is EMC
 211       fADCchanelEmc   = fCalibData->GetADCchannelEmc (module,column,row);
 212       return amp*fADCchanelEmc ;
 213     }
 214   }
 215   else{ //simulation
 216     if(absId <= fEmcCrystals) // this is EMC
 217       return fADCpedestalEmc + amp*fADCchanelEmc ;
 218   }
 219   return 0;
 220 }
 221
 222 //____________________________________________________________________________
 223 Float_t  AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
 224 {
 225   // Convert digitized CPV amplitude into charge.
 226   // Calibration parameters are taken from calibration data base for raw data,
 227   // or from digitizer parameters for simulated data.
 228
 229   if(fCalibData){
 230     Int_t relId[4];
 231     AliPHOSGetter *gime = AliPHOSGetter::Instance();
 232     gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
 233     Int_t   module = relId[0];
 234     Int_t   column = relId[3];
 235     Int_t   row    = relId[2];
 236     if(absId > fEmcCrystals) { // this is CPV
 237       fADCchanelCpv   = fCalibData->GetADCchannelCpv (module,column,row);
 238       fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
 239       return fADCpedestalCpv + amp*fADCchanelCpv ;
 240     }
 241   }
 242   else{ //simulation
 243     if(absId > fEmcCrystals) // this is CPV
 244       return fADCpedestalCpv+ amp*fADCchanelCpv ;
 245   }
 246   return 0;
 247 }
 248
 249 //____________________________________________________________________________
 250 void AliPHOSClusterizerv1::Exec(Option_t *option)
 251 {
 252   // Steering method to perform clusterization for events
 253   // in the range from fFirstEvent to fLastEvent.
 254   // This range is optionally set by SetEventRange().
 255   // if fLastEvent=-1 (by default), then process events until the end.
 256
 257   if(strstr(option,"tim"))
 258     gBenchmark->Start("PHOSClusterizer");
 259
 260   if(strstr(option,"print")) {
 261     Print() ;
 262     return ;
 263   }
 264
 265   GetCalibrationParameters() ;
 266
 267   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 268   if (fRawReader == 0)
 269     gime->SetRawDigits(kFALSE);
 270   else
 271     gime->SetRawDigits(kTRUE);
 272
 273   if (fLastEvent == -1)
 274     fLastEvent = gime->MaxEvent() - 1 ;
 275   else
 276     fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
 277   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 278
 279   Int_t ievent ;
 280   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 281     if (fRawReader == 0)
 282       gime->Event(ievent    ,"D"); // Read digits from simulated data
 283     else {
 284       gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
 285     }
 286     fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
 287
 288     MakeClusters() ;
 289
 290     AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
 291            gime->EmcRecPoints()->GetEntries(),ievent));
 292     if(AliLog::GetGlobalDebugLevel()>1)
 293       gime->EmcRecPoints()->Print();
 294
 295     if(fToUnfold)
 296       MakeUnfolding() ;
 297
 298     WriteRecPoints();
 299
 300     if(strstr(option,"deb"))
 301       PrintRecPoints(option) ;
 302
 303     //increment the total number of recpoints per run
 304     fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
 305     fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
 306   }
 307
 308   if(fWrite) //do not unload in "on flight" mode
 309     Unload();
 310
 311   if(strstr(option,"tim")){
 312     gBenchmark->Stop("PHOSClusterizer");
 313     AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
 314          gBenchmark->GetCpuTime("PHOSClusterizer"),
 315          gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
 316   }
 317 }
 318
 319 //____________________________________________________________________________
 320 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 321                                     Int_t nPar, Float_t * fitparameters) const
 322 {
 323   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 324   // The initial values for fitting procedure are set equal to the positions of local maxima.
 325   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 326
 327
 328   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 329   TClonesArray * digits = gime->Digits();
 330
 331
 332   gMinuit->mncler();                     // Reset Minuit's list of paramters
 333   gMinuit->SetPrintLevel(-1) ;           // No Printout
 334   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 335                                          // To set the address of the minimization function
 336
 337   TList * toMinuit = new TList();
 338   toMinuit->AddAt(emcRP,0) ;
 339   toMinuit->AddAt(digits,1) ;
 340
 341   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 342
 343   // filling initial values for fit parameters
 344   AliPHOSDigit * digit ;
 345
 346   Int_t ierflg  = 0;
 347   Int_t index   = 0 ;
 348   Int_t nDigits = (Int_t) nPar / 3 ;
 349
 350   Int_t iDigit ;
 351
 352   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 353
 354   for(iDigit = 0; iDigit < nDigits; iDigit++){
 355     digit = maxAt[iDigit];
 356
 357     Int_t relid[4] ;
 358     Float_t x = 0.;
 359     Float_t z = 0.;
 360     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 361     geom->RelPosInModule(relid, x, z) ;
 362
 363     Float_t energy = maxAtEnergy[iDigit] ;
 364
 365     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 366     index++ ;
 367     if(ierflg != 0){
 368       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 369       return kFALSE;
 370     }
 371     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 372     index++ ;
 373     if(ierflg != 0){
 374        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 375       return kFALSE;
 376     }
 377     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 378     index++ ;
 379     if(ierflg != 0){
 380       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 381       return kFALSE;
 382     }
 383   }
 384
 385   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 386                       //  depends on it.
 387   Double_t p1 = 1.0 ;
 388   Double_t p2 = 0.0 ;
 389
 390   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 391   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 392   gMinuit->SetMaxIterations(5);
 393   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 394
 395   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 396
 397   if(ierflg == 4){  // Minimum not found
 398     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 399     return kFALSE ;
 400   }
 401   for(index = 0; index < nPar; index++){
 402     Double_t err ;
 403     Double_t val ;
 404     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 405     fitparameters[index] = val ;
 406    }
 407
 408   delete toMinuit ;
 409   return kTRUE;
 410
 411 }
 412
 413 //____________________________________________________________________________
 414 void AliPHOSClusterizerv1::GetCalibrationParameters()
 415 {
 416   // Set calibration parameters:
 417   // if calibration database exists, they are read from database,
 418   // otherwise, they are taken from digitizer.
 419   //
 420   // It is a user responsilibity to open CDB before reconstruction, for example:
 421   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 422
 423   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 424   // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
 425
 426   fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 427
 428   if(!fCalibData)
 429     {
 430       if ( !gime->Digitizer() )
 431         gime->LoadDigitizer();
 432       AliPHOSDigitizer * dig = gime->Digitizer();
 433       fADCchanelEmc   = dig->GetEMCchannel() ;
 434       fADCpedestalEmc = dig->GetEMCpedestal();
 435
 436       fADCchanelCpv   = dig->GetCPVchannel() ;
 437       fADCpedestalCpv = dig->GetCPVpedestal() ;
 438     }
 439 }
 440
 441 //____________________________________________________________________________
 442 void AliPHOSClusterizerv1::Init()
 443 {
 444   // Make all memory allocations which can not be done in default constructor.
 445   // Attach the Clusterizer task to the list of PHOS tasks
 446
 447   AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
 448   if(!gime)
 449     gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
 450
 451   AliPHOSGeometry * geom = gime->PHOSGeometry();
 452
 453   fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
 454
 455   if(!gMinuit)
 456     gMinuit = new TMinuit(100);
 457
 458   if ( !gime->Clusterizer() ) {
 459     gime->PostClusterizer(this);
 460   }
 461 }
 462
 463 //____________________________________________________________________________
 464 void AliPHOSClusterizerv1::InitParameters()
 465 {
 466
 467   fNumberOfCpvClusters     = 0 ;
 468   fNumberOfEmcClusters     = 0 ;
 469
 470   fCpvClusteringThreshold  = 0.0;
 471   fEmcClusteringThreshold  = 0.2;
 472
 473   fEmcLocMaxCut            = 0.03 ;
 474   fCpvLocMaxCut            = 0.03 ;
 475
 476   fEmcMinE                 = 0.01 ;
 477   fCpvMinE                 = 0.0  ;
 478
 479   fW0                      = 4.5 ;
 480   fW0CPV                   = 4.0 ;
 481
 482   fEmcTimeGate             = 1.e-8 ;
 483
 484   fToUnfold                = kTRUE ;
 485
 486   fRecPointsInRun          = 0 ;
 487
 488   fWrite                   = kTRUE ;
 489
 490   fCalibData               = 0 ;
 491
 492   SetEventRange(0,-1) ;
 493
 494   fIsOldRCUFormat          = kFALSE;
 495 }
 496
 497 //____________________________________________________________________________
 498 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 499 {
 500   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 501   //                                       = 1 are neighbour
 502   //                                       = 2 are not neighbour but do not continue searching
 503   // neighbours are defined as digits having at least a common vertex
 504   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 505   //                                      which is compared to a digit (d2)  not yet in a cluster
 506
 507   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 508
 509   Int_t rv = 0 ;
 510
 511   Int_t relid1[4] ;
 512   geom->AbsToRelNumbering(d1->GetId(), relid1) ;
 513
 514   Int_t relid2[4] ;
 515   geom->AbsToRelNumbering(d2->GetId(), relid2) ;
 516
 517   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 518     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 519     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 520
 521     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 522       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 523       rv = 1 ;
 524     }
 525     else {
 526       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 527         rv = 2; //  Difference in row numbers is too large to look further
 528     }
 529
 530   }
 531   else {
 532
 533     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 534       rv=2 ;
 535
 536   }
 537
 538   return rv ;
 539 }
 540 //____________________________________________________________________________
 541 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
 542 {
 543   // Remove digits with amplitudes below threshold
 544
 545   for(Int_t i=0; i<digits->GetEntriesFast(); i++){
 546     AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
 547     if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
 548          (IsInCpv(digit) && CalibrateCPV(digit->GetAmp()   ,digit->GetId()) < fCpvMinE) )
 549       digits->RemoveAt(i) ;
 550   }
 551   digits->Compress() ;
 552   for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
 553     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 554     digit->SetIndexInList(i) ;
 555   }
 556
 557 }
 558 //____________________________________________________________________________
 559 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 560 {
 561   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 562
 563   Bool_t rv = kFALSE ;
 564   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 565
 566   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 567
 568   if(digit->GetId() <= nEMC )   rv = kTRUE;
 569
 570   return rv ;
 571 }
 572
 573 //____________________________________________________________________________
 574 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 575 {
 576   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 577
 578   Bool_t rv = kFALSE ;
 579
 580   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 581
 582   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 583
 584   if(digit->GetId() > nEMC )   rv = kTRUE;
 585
 586   return rv ;
 587 }
 588
 589 //____________________________________________________________________________
 590 void AliPHOSClusterizerv1::Unload()
 591 {
 592   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 593   gime->PhosLoader()->UnloadDigits() ;
 594   gime->PhosLoader()->UnloadRecPoints() ;
 595 }
 596
 597 //____________________________________________________________________________
 598 void AliPHOSClusterizerv1::WriteRecPoints()
 599 {
 600
 601   // Creates new branches with given title
 602   // fills and writes into TreeR.
 603
 604   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 605
 606   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 607   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 608   TClonesArray * digits = gime->Digits() ;
 609
 610
 611   Int_t index ;
 612   //Evaluate position, dispersion and other RecPoint properties..
 613   Int_t nEmc = emcRecPoints->GetEntriesFast();
 614   for(index = 0; index < nEmc; index++){
 615     AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
 616     rp->Purify(fEmcMinE) ;
 617     if(rp->GetMultiplicity()>0.) //If this RP is not empty
 618       rp->EvalAll(fW0,digits) ;
 619     else{
 620       emcRecPoints->RemoveAt(index) ;
 621       delete rp ;
 622     }
 623   }
 624   emcRecPoints->Compress() ;
 625   emcRecPoints->Sort() ;
 626   //  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
 627   for(index = 0; index < emcRecPoints->GetEntries(); index++){
 628     dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
 629   }
 630
 631   //Now the same for CPV
 632   for(index = 0; index < cpvRecPoints->GetEntries(); index++){
 633     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
 634     rp->EvalAll(fW0CPV,digits) ;
 635   }
 636   cpvRecPoints->Sort() ;
 637
 638   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 639     dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
 640
 641   cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
 642
 643   if(fWrite){ //We write TreeR
 644     TTree * treeR = gime->TreeR();
 645
 646     Int_t bufferSize = 32000 ;
 647     Int_t splitlevel = 0 ;
 648
 649     //First EMC
 650     TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
 651     emcBranch->SetTitle(BranchName());
 652
 653     //Now CPV branch
 654     TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
 655     cpvBranch->SetTitle(BranchName());
 656
 657     emcBranch        ->Fill() ;
 658     cpvBranch        ->Fill() ;
 659
 660     gime->WriteRecPoints("OVERWRITE");
 661     gime->WriteClusterizer("OVERWRITE");
 662   }
 663 }
 664
 665 //____________________________________________________________________________
 666 void AliPHOSClusterizerv1::MakeClusters()
 667 {
 668   // Steering method to construct the clusters stored in a list of Reconstructed Points
 669   // A cluster is defined as a list of neighbour digits
 670
 671
 672   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 673
 674   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 675   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 676
 677   emcRecPoints->Delete() ;
 678   cpvRecPoints->Delete() ;
 679
 680   TClonesArray * digits = gime->Digits() ;
 681
 682   //Remove digits below threshold
 683   CleanDigits(digits) ;
 684
 685   TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
 686
 687   // Clusterization starts
 688
 689   TIter nextdigit(digitsC) ;
 690   AliPHOSDigit * digit ;
 691   Bool_t notremoved = kTRUE ;
 692
 693   while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
 694
 695
 696     AliPHOSRecPoint * clu = 0 ;
 697
 698     TArrayI clusterdigitslist(1500) ;
 699     Int_t index ;
 700
 701     if (( IsInEmc (digit) &&
 702           CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
 703         ( IsInCpv (digit) &&
 704           CalibrateCPV(digit->GetAmp()   ,digit->GetId()) > fCpvClusteringThreshold ) ) {
 705       Int_t iDigitInCluster = 0 ;
 706
 707       if  ( IsInEmc(digit) ) {
 708         // start a new EMC RecPoint
 709         if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 710           emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 711
 712         emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 713         clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
 714         fNumberOfEmcClusters++ ;
 715         clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
 716         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 717         iDigitInCluster++ ;
 718         digitsC->Remove(digit) ;
 719
 720       } else {
 721
 722         // start a new CPV cluster
 723         if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 724           cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
 725
 726         cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 727
 728         clu =  dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
 729         fNumberOfCpvClusters++ ;
 730         clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
 731         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 732         iDigitInCluster++ ;
 733         digitsC->Remove(digit) ;
 734         nextdigit.Reset() ;
 735
 736         // Here we remove remaining EMC digits, which cannot make a cluster
 737
 738         if( notremoved ) {
 739           while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
 740             if( IsInEmc(digit) )
 741               digitsC->Remove(digit) ;
 742             else
 743               break ;
 744           }
 745           notremoved = kFALSE ;
 746         }
 747
 748       } // else
 749
 750       nextdigit.Reset() ;
 751
 752       AliPHOSDigit * digitN ;
 753       index = 0 ;
 754       while (index < iDigitInCluster){ // scan over digits already in cluster
 755         digit =  dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) )  ;
 756         index++ ;
 757         while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
 758           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 759           switch (ineb ) {
 760           case 0 :   // not a neighbour
 761             break ;
 762           case 1 :   // are neighbours
 763             if (IsInEmc (digitN))
 764               clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
 765             else
 766               clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp()   , digitN->GetId() ) );
 767
 768             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
 769             iDigitInCluster++ ;
 770             digitsC->Remove(digitN) ;
 771             break ;
 772           case 2 :   // too far from each other
 773             goto endofloop;
 774           } // switch
 775
 776         } // while digitN
 777
 778       endofloop: ;
 779         nextdigit.Reset() ;
 780
 781       } // loop over cluster
 782
 783     } // energy theshold
 784
 785
 786   } // while digit
 787
 788   delete digitsC ;
 789
 790 }
 791
 792 //____________________________________________________________________________
 793 void AliPHOSClusterizerv1::MakeUnfolding()
 794 {
 795   // Unfolds clusters using the shape of an ElectroMagnetic shower
 796   // Performs unfolding of all EMC/CPV clusters
 797
 798   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 799
 800   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 801
 802   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 803   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 804   TClonesArray * digits = gime->Digits() ;
 805
 806   // Unfold first EMC clusters
 807   if(fNumberOfEmcClusters > 0){
 808
 809     Int_t nModulesToUnfold = geom->GetNModules() ;
 810
 811     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 812     Int_t index ;
 813     for(index = 0 ; index < numberofNotUnfolded ; index++){
 814
 815       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
 816       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 817         break ;
 818
 819       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 820       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 821       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 822       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
 823
 824       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 825         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 826         emcRecPoints->Remove(emcRecPoint);
 827         emcRecPoints->Compress() ;
 828         index-- ;
 829         fNumberOfEmcClusters -- ;
 830         numberofNotUnfolded-- ;
 831       }
 832       else{
 833         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 834       }
 835
 836       delete[] maxAt ;
 837       delete[] maxAtEnergy ;
 838     }
 839   }
 840   // Unfolding of EMC clusters finished
 841
 842
 843   // Unfold now CPV clusters
 844   if(fNumberOfCpvClusters > 0){
 845
 846     Int_t nModulesToUnfold = geom->GetNModules() ;
 847
 848     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 849     Int_t index ;
 850     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 851
 852       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
 853
 854       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 855         break ;
 856
 857       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 858
 859       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 860       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 861       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 862       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
 863
 864       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 865         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 866         cpvRecPoints->Remove(emcRecPoint);
 867         cpvRecPoints->Compress() ;
 868         index-- ;
 869         numberofCpvNotUnfolded-- ;
 870         fNumberOfCpvClusters-- ;
 871       }
 872
 873       delete[] maxAt ;
 874       delete[] maxAtEnergy ;
 875     }
 876   }
 877   //Unfolding of Cpv clusters finished
 878
 879 }
 880
 881 //____________________________________________________________________________
 882 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 883 {
 884   // Shape of the shower (see PHOS TDR)
 885   // If you change this function, change also the gradient evaluation in ChiSquare()
 886
 887   Double_t r4    = r*r*r*r ;
 888   Double_t r295  = TMath::Power(r, 2.95) ;
 889   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 890   return shape ;
 891 }
 892
 893 //____________________________________________________________________________
 894 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 895                                           Int_t nMax,
 896                                           AliPHOSDigit ** maxAt,
 897                                           Float_t * maxAtEnergy)
 898 {
 899   // Performs the unfolding of a cluster with nMax overlapping showers
 900
 901   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 902
 903   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 904
 905   const TClonesArray * digits = gime->Digits() ;
 906   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 907   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 908
 909   Int_t nPar = 3 * nMax ;
 910   Float_t * fitparameters = new Float_t[nPar] ;
 911
 912   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 913   if( !rv ) {
 914     // Fit failed, return and remove cluster
 915     iniEmc->SetNExMax(-1) ;
 916     delete[] fitparameters ;
 917     return ;
 918   }
 919
 920   // create ufolded rec points and fill them with new energy lists
 921   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 922   // and later correct this number in acordance with actual energy deposition
 923
 924   Int_t nDigits = iniEmc->GetMultiplicity() ;
 925   Float_t * efit = new Float_t[nDigits] ;
 926   Float_t xDigit=0.,zDigit=0.,distance=0. ;
 927   Float_t xpar=0.,zpar=0.,epar=0.  ;
 928   Int_t relid[4] ;
 929   AliPHOSDigit * digit = 0 ;
 930   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 931
 932   Int_t iparam ;
 933   Int_t iDigit ;
 934   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 935     digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
 936     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 937     geom->RelPosInModule(relid, xDigit, zDigit) ;
 938     efit[iDigit] = 0;
 939
 940     iparam = 0 ;
 941     while(iparam < nPar ){
 942       xpar = fitparameters[iparam] ;
 943       zpar = fitparameters[iparam+1] ;
 944       epar = fitparameters[iparam+2] ;
 945       iparam += 3 ;
 946       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 947       distance =  TMath::Sqrt(distance) ;
 948       efit[iDigit] += epar * ShowerShape(distance) ;
 949     }
 950   }
 951
 952
 953   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 954   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 955   // to its contribution to efit
 956
 957   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 958   Float_t ratio ;
 959
 960   iparam = 0 ;
 961   while(iparam < nPar ){
 962     xpar = fitparameters[iparam] ;
 963     zpar = fitparameters[iparam+1] ;
 964     epar = fitparameters[iparam+2] ;
 965     iparam += 3 ;
 966
 967     AliPHOSEmcRecPoint * emcRP = 0 ;
 968
 969     if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
 970
 971       if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 972         emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
 973
 974       (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 975       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
 976       fNumberOfEmcClusters++ ;
 977       emcRP->SetNExMax((Int_t)nPar/3) ;
 978     }
 979     else{//create new entries in fCpvRecPoints
 980       if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 981         cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
 982
 983       (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 984       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
 985       fNumberOfCpvClusters++ ;
 986     }
 987
 988     Float_t eDigit ;
 989     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 990       digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
 991       geom->AbsToRelNumbering(digit->GetId(), relid) ;
 992       geom->RelPosInModule(relid, xDigit, zDigit) ;
 993       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 994       distance =  TMath::Sqrt(distance) ;
 995       ratio = epar * ShowerShape(distance) / efit[iDigit] ;
 996       eDigit = emcEnergies[iDigit] * ratio ;
 997       emcRP->AddDigit( *digit, eDigit ) ;
 998     }
 999   }
1000
1001   delete[] fitparameters ;
1002   delete[] efit ;
1003
1004 }
1005
1006 //_____________________________________________________________________________
1007 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1008 {
1009   // Calculates the Chi square for the cluster unfolding minimization
1010   // Number of parameters, Gradient, Chi squared, parameters, what to do
1011
1012   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
1013
1014   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
1015   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
1016
1017
1018
1019   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
1020
1021   Int_t * emcDigits     = emcRP->GetDigitsList() ;
1022
1023   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1024
1025   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1026
1027   const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1028   fret = 0. ;
1029   Int_t iparam ;
1030
1031   if(iflag == 2)
1032     for(iparam = 0 ; iparam < nPar ; iparam++)
1033       Grad[iparam] = 0 ; // Will evaluate gradient
1034
1035   Double_t efit ;
1036
1037   AliPHOSDigit * digit ;
1038   Int_t iDigit ;
1039
1040   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1041
1042     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1043
1044     Int_t relid[4] ;
1045     Float_t xDigit ;
1046     Float_t zDigit ;
1047
1048     geom->AbsToRelNumbering(digit->GetId(), relid) ;
1049
1050     geom->RelPosInModule(relid, xDigit, zDigit) ;
1051
1052      if(iflag == 2){  // calculate gradient
1053        Int_t iParam = 0 ;
1054        efit = 0 ;
1055        while(iParam < nPar ){
1056          Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
1057          iParam++ ;
1058          distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
1059          distance = TMath::Sqrt( distance ) ;
1060          iParam++ ;
1061          efit += x[iParam] * ShowerShape(distance) ;
1062          iParam++ ;
1063        }
1064        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1065        iParam = 0 ;
1066        while(iParam < nPar ){
1067          Double_t xpar = x[iParam] ;
1068          Double_t zpar = x[iParam+1] ;
1069          Double_t epar = x[iParam+2] ;
1070          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1071          Double_t shape = sum * ShowerShape(dr) ;
1072          Double_t r4 = dr*dr*dr*dr ;
1073          Double_t r295 = TMath::Power(dr,2.95) ;
1074          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1075                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1076
1077          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
1078          iParam++ ;
1079          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
1080          iParam++ ;
1081          Grad[iParam] += shape ;                                  // Derivative over energy
1082          iParam++ ;
1083        }
1084      }
1085      efit = 0;
1086      iparam = 0 ;
1087
1088      while(iparam < nPar ){
1089        Double_t xpar = x[iparam] ;
1090        Double_t zpar = x[iparam+1] ;
1091        Double_t epar = x[iparam+2] ;
1092        iparam += 3 ;
1093        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1094        distance =  TMath::Sqrt(distance) ;
1095        efit += epar * ShowerShape(distance) ;
1096      }
1097
1098      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1099      // Here we assume, that sigma = sqrt(E)
1100   }
1101
1102 }
1103
1104 //____________________________________________________________________________
1105 void AliPHOSClusterizerv1::Print(const Option_t *)const
1106 {
1107   // Print clusterizer parameters
1108
1109   TString message ;
1110   TString taskName(GetName()) ;
1111   taskName.ReplaceAll(Version(), "") ;
1112
1113   if( strcmp(GetName(), "") !=0 ) {
1114     // Print parameters
1115     message  = "\n--------------- %s %s -----------\n" ;
1116     message += "Clusterizing digits from the file: %s\n" ;
1117     message += "                           Branch: %s\n" ;
1118     message += "                       EMC Clustering threshold = %f\n" ;
1119     message += "                       EMC Local Maximum cut    = %f\n" ;
1120     message += "                       EMC Logarothmic weight   = %f\n" ;
1121     message += "                       CPV Clustering threshold = %f\n" ;
1122     message += "                       CPV Local Maximum cut    = %f\n" ;
1123     message += "                       CPV Logarothmic weight   = %f\n" ;
1124     if(fToUnfold)
1125       message += " Unfolding on\n" ;
1126     else
1127       message += " Unfolding off\n" ;
1128
1129     message += "------------------------------------------------------------------" ;
1130   }
1131   else
1132     message = " AliPHOSClusterizerv1 not initialized " ;
1133
1134   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1135        taskName.Data(),
1136        GetTitle(),
1137        taskName.Data(),
1138        GetName(),
1139        fEmcClusteringThreshold,
1140        fEmcLocMaxCut,
1141        fW0,
1142        fCpvClusteringThreshold,
1143        fCpvLocMaxCut,
1144        fW0CPV )) ;
1145 }
1146
1147
1148 //____________________________________________________________________________
1149 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1150 {
1151   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1152
1153   AliPHOSGetter * gime = AliPHOSGetter::Instance();
1154
1155   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1156   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1157
1158   AliInfo(Form("\nevent %d \n    Found %d EMC RecPoints and %d CPV RecPoints",
1159                gAlice->GetEvNumber(),
1160                emcRecPoints->GetEntriesFast(),
1161                cpvRecPoints->GetEntriesFast() ))  ;
1162
1163   fRecPointsInRun +=  emcRecPoints->GetEntriesFast() ;
1164   fRecPointsInRun +=  cpvRecPoints->GetEntriesFast() ;
1165
1166
1167   if(strstr(option,"all")) {
1168     printf("\n EMC clusters \n") ;
1169     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1170     Int_t index ;
1171     for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1172       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1173       TVector3  locpos;
1174       rp->GetLocalPosition(locpos);
1175       Float_t lambda[2];
1176       rp->GetElipsAxis(lambda);
1177       Int_t * primaries;
1178       Int_t nprimaries;
1179       primaries = rp->GetPrimaries(nprimaries);
1180       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1181               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1182               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1183
1184       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1185         printf("%d ", primaries[iprimary] ) ;
1186       }
1187       printf("\n") ;
1188     }
1189
1190     //Now plot CPV recPoints
1191     printf("\n CPV clusters \n") ;
1192     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1193     for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1194       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1195
1196       TVector3  locpos;
1197       rp->GetLocalPosition(locpos);
1198
1199       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1200              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1201              locpos.X(), locpos.Y(), locpos.Z()) ;
1202     }
1203   }
1204 }
1205