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