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