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