Introducing Riostream.h
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18
19 /* $Log:
20
21    1 October 2000. Yuri Kharlov:
22
23      AreNeighbours()
24
25      PPSD upper layer is considered if number of layers>1
26
27
28
29    18 October 2000. Yuri Kharlov:
30
31      AliEMCALClusterizerv1()
32
33      CPV clusterizing parameters added
34
35
36
37      MakeClusters()
38
39      After first PPSD digit remove EMC digits only once
40
41 */
42
43 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
44
45 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
46
47 //                           of new  IO (à la PHOS)
48
49 //////////////////////////////////////////////////////////////////////////////
50
51 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
52
53 //  unfolds the clusters having several local maxima.  
54
55 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
56
57 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all 
58
59 //  parameters including input digits branch title, thresholds etc.)
60
61 //  This TTask is normally called from Reconstructioner, but can as well be used in 
62
63 //  standalone mode.
64
65 // Use Case:
66
67 //  root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")  
68
69 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
70
71 //               //reads gAlice from header file "..."                      
72
73 //  root [1] cl->ExecuteTask()  
74
75 //               //finds RecPoints in all events stored in galice.root
76
77 //  root [2] cl->SetDigitsBranch("digits2") 
78
79 //               //sets another title for Digitis (input) branch
80
81 //  root [3] cl->SetRecPointsBranch("recp2")  
82
83 //               //sets another title four output branches
84
85 //  root [4] cl->SetTowerLocalMaxCut(0.03)  
86
87 //               //set clusterization parameters
88
89 //  root [5] cl->ExecuteTask("deb all time")  
90
91 //               //once more finds RecPoints options are 
92
93 //               // deb - print number of found rec points
94
95 //               // deb all - print number of found RecPoints and some their characteristics 
96
97 //               // time - print benchmarking results
98
99
100
101 // --- ROOT system ---
102
103
104
105 #include "TROOT.h" 
106
107 #include "TFile.h" 
108
109 #include "TFolder.h" 
110
111 #include "TMath.h" 
112
113 #include "TMinuit.h"
114
115 #include "TTree.h" 
116
117 #include "TSystem.h" 
118
119 #include "TBenchmark.h"
120
121
122
123 // --- Standard library ---
124
125
126
127 #include <Riostream.h>
128
129
130
131 // --- AliRoot header files ---
132
133
134
135 #include "AliEMCALClusterizerv1.h"
136
137 #include "AliEMCALDigit.h"
138
139 #include "AliEMCALDigitizer.h"
140
141 #include "AliEMCALTowerRecPoint.h"
142
143 #include "AliEMCAL.h"
144
145 #include "AliEMCALGetter.h"
146
147 #include "AliEMCALGeometry.h"
148
149 #include "AliRun.h"
150
151
152
153 ClassImp(AliEMCALClusterizerv1)
154
155   
156
157 //____________________________________________________________________________
158
159   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
160
161 {
162
163   // default ctor (to be used mainly by Streamer)
164
165   
166
167   InitParameters() ; 
168
169   fDefaultInit = kTRUE ; 
170
171 }
172
173
174
175 //____________________________________________________________________________
176
177 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
178
179 :AliEMCALClusterizer(headerFile, name, toSplit)
180
181 {
182
183   // ctor with the indication of the file where header Tree and digits Tree are stored
184
185   
186
187   InitParameters() ; 
188
189   Init() ;
190
191   fDefaultInit = kFALSE ; 
192
193
194
195 }
196
197
198
199 //____________________________________________________________________________
200
201   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
202
203 {
204
205   // dtor
206
207   fSplitFile = 0 ; 
208
209   
210
211 }
212
213
214
215 //____________________________________________________________________________
216
217 const TString AliEMCALClusterizerv1::BranchName() const 
218
219 {  
220
221   TString branchName(GetName() ) ;
222
223   branchName.Remove(branchName.Index(Version())-1) ;
224
225   return branchName ;
226
227 }
228
229
230
231 //____________________________________________________________________________
232
233 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const
234
235 {//To be replased later by the method, reading individual parameters from the database
236
237
238
239   if ( inpresho ) // calibrate as pre shower
240
241      return -fADCpedestalPreSho + amp * fADCchannelPreSho ; 
242
243
244
245   else //calibrate as tower 
246
247     return -fADCpedestalTower + amp * fADCchannelTower ;                
248
249 }
250
251
252
253 //____________________________________________________________________________
254
255 void AliEMCALClusterizerv1::Exec(Option_t * option)
256
257 {
258
259   // Steering method
260
261
262
263   if( strcmp(GetName(), "")== 0 ) 
264
265     Init() ;
266
267
268
269   if(strstr(option,"tim"))
270
271     gBenchmark->Start("EMCALClusterizer"); 
272
273   
274
275   if(strstr(option,"print"))
276
277     Print("") ; 
278
279
280
281   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
282
283   if(gime->BranchExists("RecPoints"))
284
285     return ;
286
287   Int_t nevents = gime->MaxEvent() ;
288
289   Int_t ievent ;
290
291
292
293   for(ievent = 0; ievent < nevents; ievent++){
294
295
296
297     gime->Event(ievent,"D") ;
298
299
300
301     if(ievent == 0)
302
303       GetCalibrationParameters() ;
304
305
306
307     fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ;
308
309            
310
311     MakeClusters() ;
312
313     
314
315     if(fToUnfold)
316
317       MakeUnfolding() ;
318
319
320
321     WriteRecPoints(ievent) ;
322
323
324
325     if(strstr(option,"deb"))  
326
327       PrintRecPoints(option) ;
328
329
330
331     //increment the total number of digits per run 
332
333     fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ;  
334
335     fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ;  
336
337  }
338
339   
340
341   if(strstr(option,"tim")){
342
343     gBenchmark->Stop("EMCALClusterizer");
344
345     cout << "AliEMCALClusterizer:" << endl ;
346
347     cout << "  took " << gBenchmark->GetCpuTime("EMCALClusterizer") << " seconds for Clusterizing " 
348
349          <<  gBenchmark->GetCpuTime("EMCALClusterizer")/nevents << " seconds per event " << endl ;
350
351     cout << endl ;
352
353   }
354
355   
356
357 }
358
359
360
361 //____________________________________________________________________________
362
363 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
364
365                                     Int_t nPar, Float_t * fitparameters) const
366
367
368
369   // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
370
371   // The initial values for fitting procedure are set equal to the positions of local maxima.
372
373   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
374
375
376
377   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
378
379   TClonesArray * digits = gime->Digits() ; 
380
381   
382
383
384
385   gMinuit->mncler();                     // Reset Minuit's list of paramters
386
387   gMinuit->SetPrintLevel(-1) ;           // No Printout
388
389   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;  
390
391                                          // To set the address of the minimization function 
392
393
394
395   TList * toMinuit = new TList();
396
397   toMinuit->AddAt(emcRP,0) ;
398
399   toMinuit->AddAt(digits,1) ;
400
401   
402
403   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
404
405
406
407   // filling initial values for fit parameters
408
409   AliEMCALDigit * digit ;
410
411
412
413   Int_t ierflg  = 0; 
414
415   Int_t index   = 0 ;
416
417   Int_t nDigits = (Int_t) nPar / 3 ;
418
419
420
421   Int_t iDigit ;
422
423
424
425   AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
426
427
428
429   for(iDigit = 0; iDigit < nDigits; iDigit++){
430
431     digit = maxAt[iDigit]; 
432
433
434
435     Int_t relid[4] ;
436
437     Float_t x = 0.;
438
439     Float_t z = 0.;
440
441     geom->AbsToRelNumbering(digit->GetId(), relid) ;
442
443     geom->PosInAlice(relid, x, z) ;
444
445
446
447     Float_t energy = maxAtEnergy[iDigit] ;
448
449
450
451     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
452
453     index++ ;   
454
455     if(ierflg != 0){ 
456
457       cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
458
459       return kFALSE;
460
461     }
462
463     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
464
465     index++ ;   
466
467     if(ierflg != 0){
468
469       cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
470
471       return kFALSE;
472
473     }
474
475     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
476
477     index++ ;   
478
479     if(ierflg != 0){
480
481       cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : energy = " << energy << endl ;      
482
483       return kFALSE;
484
485     }
486
487   }
488
489
490
491   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
492
493                       //  depends on it. 
494
495   Double_t p1 = 1.0 ;
496
497   Double_t p2 = 0.0 ;
498
499
500
501   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
502
503   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
504
505   gMinuit->SetMaxIterations(5);
506
507   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
508
509
510
511   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
512
513
514
515   if(ierflg == 4){  // Minimum not found   
516
517     cout << "EMCAL Unfolding>  Fit not converged, cluster abandoned "<< endl ;      
518
519     return kFALSE ;
520
521   }            
522
523   for(index = 0; index < nPar; index++){
524
525     Double_t err ;
526
527     Double_t val ;
528
529     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
530
531     fitparameters[index] = val ;
532
533    }
534
535
536
537   delete toMinuit ;
538
539   return kTRUE;
540
541
542
543 }
544
545
546
547 //____________________________________________________________________________
548
549 void AliEMCALClusterizerv1::GetCalibrationParameters() 
550
551 {
552
553   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
554
555   const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
556
557
558
559   fADCchannelTower   = dig->GetTowerchannel() ;
560
561   fADCpedestalTower  = dig->GetTowerpedestal();
562
563
564
565   fADCchannelPreSho  = dig->GetPreShochannel() ;
566
567   fADCpedestalPreSho = dig->GetPreShopedestal() ; 
568
569
570
571 }
572
573
574
575 //____________________________________________________________________________
576
577 void AliEMCALClusterizerv1::Init()
578
579 {
580
581   // Make all memory allocations which can not be done in default constructor.
582
583   // Attach the Clusterizer task to the list of EMCAL tasks
584
585   
586
587   if ( strcmp(GetTitle(), "") == 0 )
588
589     SetTitle("galice.root") ;
590
591
592
593   TString branchname = GetName() ;
594
595   branchname.Remove(branchname.Index(Version())-1) ;
596
597
598
599   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ; 
600
601   if ( gime == 0 ) {
602
603     cerr << "ERROR: AliEMCALClusterizerv1::Init -> Could not obtain the Getter object !" << endl ; 
604
605     return ;
606
607   } 
608
609
610
611   fSplitFile = 0 ;
612
613   if(fToSplit){
614
615     // construct the name of the file as /path/EMCAL.SDigits.root
616
617     //First - extract full path if necessary
618
619     TString fileName(GetTitle()) ;
620
621     Ssiz_t islash = fileName.Last('/') ;
622
623     if(islash<fileName.Length())
624
625       fileName.Remove(islash+1,fileName.Length()) ;
626
627     else
628
629       fileName="" ;
630
631     // Next - append the file name 
632
633     fileName+="EMCAL.RecData." ;
634
635     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
636
637       fileName+=branchname ;
638
639       fileName+="." ;
640
641     }
642
643     fileName+="root" ;
644
645     // Finally - check if the file already opened or open the file
646
647     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));   
648
649     if(!fSplitFile)
650
651       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
652
653   }
654
655
656
657     
658
659   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
660
661   fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
662
663
664
665   if(!gMinuit) 
666
667     gMinuit = new TMinuit(100) ;
668
669
670
671   gime->PostClusterizer(this) ;
672
673   gime->PostRecPoints(branchname ) ;
674
675  
676
677 }
678
679
680
681 //____________________________________________________________________________
682
683 void AliEMCALClusterizerv1::InitParameters()
684
685 {
686
687   fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ; 
688
689
690
691  
692
693   
694
695   fPreShoClusteringThreshold  = 0.0001;
696
697   fTowerClusteringThreshold   = 0.2;   
698
699   
700
701   fTowerLocMaxCut  = 0.03 ;
702
703   fPreShoLocMaxCut = 0.03 ;
704
705   
706
707   fW0     = 4.5 ;
708
709   fW0CPV  = 4.0 ;
710
711
712
713   fTimeGate = 1.e-8 ; 
714
715   
716
717   fToUnfold = kFALSE ;
718
719    
720
721   TString clusterizerName( GetName()) ; 
722
723   if (clusterizerName.IsNull() ) 
724
725     clusterizerName = "Default" ; 
726
727   clusterizerName.Append(":") ; 
728
729   clusterizerName.Append(Version()) ; 
730
731   SetName(clusterizerName) ;
732
733   fRecPointsInRun          = 0 ; 
734
735
736
737 }
738
739
740
741 //____________________________________________________________________________
742
743 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
744
745 {
746
747   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
748
749   //                                       = 1 are neighbour
750
751   //                                       = 2 are not neighbour but do not continue searching
752
753   // neighbours are defined as digits having at least a common vertex 
754
755   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
756
757   //                                      which is compared to a digit (d2)  not yet in a cluster  
758
759
760
761    AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
762
763
764
765   Int_t rv = 0 ; 
766
767
768
769   Int_t relid1[4] ; 
770
771   geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
772
773
774
775   Int_t relid2[4] ; 
776
777   geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
778
779  
780
781   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm 
782
783     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
784
785     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
786
787     
788
789     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
790
791       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
792
793       rv = 1 ; 
794
795     }
796
797     else {
798
799       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
800
801         rv = 2; //  Difference in row numbers is too large to look further 
802
803     }
804
805
806
807   } 
808
809   else {
810
811     
812
813     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
814
815       rv=2 ;
816
817
818
819   }
820
821
822
823   return rv ; 
824
825 }
826
827
828
829
830
831 //____________________________________________________________________________
832
833 Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const
834
835 {
836
837   // Tells if (true) or not (false) the digit is in a EMCAL-Tower 
838
839  
840
841   Bool_t rv = kFALSE ; 
842
843   if (!digit->IsInPreShower()) 
844
845     rv = kTRUE; 
846
847   return rv ; 
848
849 }
850
851
852
853 //____________________________________________________________________________
854
855 Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const
856
857 {
858
859   // Tells if (true) or not (false) the digit is in a EMCAL-PreShower
860
861  
862
863   Bool_t rv = kFALSE ; 
864
865   if (digit->IsInPreShower()) 
866
867     rv = kTRUE; 
868
869   return rv ; 
870
871 }
872
873
874
875 //____________________________________________________________________________
876
877 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
878
879 {
880
881
882
883   // Creates new branches with given title
884
885   // fills and writes into TreeR.
886
887
888
889   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
890
891   TObjArray * towerRecPoints = gime->TowerRecPoints() ; 
892
893   TObjArray * preshoRecPoints = gime->PreShowerRecPoints() ; 
894
895   TClonesArray * digits = gime->Digits() ; 
896
897   TTree * treeR ; 
898
899   
900
901   if(fToSplit){
902
903     if(!fSplitFile)
904
905       return ;
906
907     fSplitFile->cd() ;
908
909     TString name("TreeR") ;
910
911     name += event ; 
912
913     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
914
915   }
916
917   else{
918
919     treeR = gAlice->TreeR();
920
921   }
922
923
924
925   if(!treeR){
926
927     gAlice->MakeTree("R", fSplitFile);
928
929     treeR = gAlice->TreeR() ;
930
931   }
932
933  
934
935   Int_t index ;
936
937   //Evaluate position, dispersion and other RecPoint properties...
938
939   for(index = 0; index < towerRecPoints->GetEntries(); index++)
940
941     (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->EvalAll(fW0,digits) ;
942
943   
944
945   towerRecPoints->Sort() ;
946
947
948
949   for(index = 0; index < towerRecPoints->GetEntries(); index++)
950
951     (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->SetIndexInList(index) ;
952
953
954
955   towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ; 
956
957
958
959   //Now the same for pre shower
960
961   for(index = 0; index < preshoRecPoints->GetEntries(); index++)
962
963     (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits)  ;
964
965
966
967   preshoRecPoints->Sort() ;
968
969
970
971   for(index = 0; index < preshoRecPoints->GetEntries(); index++)
972
973     (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->SetIndexInList(index) ;
974
975
976
977   preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ;
978
979   
980
981   Int_t bufferSize = 32000 ;    
982
983   Int_t splitlevel = 0 ;
984
985
986
987   //First Tower branch
988
989   TBranch * towerBranch = treeR->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel);
990
991   towerBranch->SetTitle(BranchName());
992
993   
994
995   //Now Pre Shower branch 
996
997   TBranch * preshoBranch = treeR->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel);
998
999   preshoBranch->SetTitle(BranchName());
1000
1001     
1002
1003   //And Finally  clusterizer branch
1004
1005   AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
1006
1007   TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
1008
1009                                               &cl,bufferSize,splitlevel);
1010
1011   clusterizerBranch->SetTitle(BranchName());
1012
1013
1014
1015   towerBranch        ->Fill() ;
1016
1017   preshoBranch        ->Fill() ;
1018
1019   clusterizerBranch->Fill() ;
1020
1021
1022
1023   treeR->AutoSave() ; //Write(0,kOverwrite) ;  
1024
1025   if(gAlice->TreeR()!=treeR)
1026
1027     treeR->Delete(); 
1028
1029 }
1030
1031
1032
1033 //____________________________________________________________________________
1034
1035 void AliEMCALClusterizerv1::MakeClusters()
1036
1037 {
1038
1039   // Steering method to construct the clusters stored in a list of Reconstructed Points
1040
1041   // A cluster is defined as a list of neighbour digits
1042
1043     
1044
1045   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
1046
1047   
1048
1049   TObjArray * towerRecPoints  = gime->TowerRecPoints(BranchName()) ; 
1050
1051   TObjArray * preshoRecPoints = gime->PreShowerRecPoints(BranchName()) ; 
1052
1053   towerRecPoints->Delete() ;
1054
1055   preshoRecPoints->Delete() ;
1056
1057   
1058
1059   TClonesArray * digits = gime->Digits() ; 
1060
1061   if ( !digits ) {
1062
1063     cerr << "ERROR:  AliEMCALClusterizerv1::MakeClusters -> Digits with name " 
1064
1065          << GetName() << " not found ! " << endl ; 
1066
1067     abort() ; 
1068
1069   } 
1070
1071   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
1072
1073   
1074
1075   
1076
1077   // Clusterization starts  
1078
1079   
1080
1081   TIter nextdigit(digitsC) ; 
1082
1083   AliEMCALDigit * digit ; 
1084
1085   Bool_t notremoved = kTRUE ;
1086
1087   
1088
1089   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
1090
1091     AliEMCALRecPoint * clu = 0 ; 
1092
1093     
1094
1095     TArrayI clusterdigitslist(1500) ;   
1096
1097     Int_t index ;
1098
1099  
1100
1101     if (( IsInTower (digit)  && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold  ) || 
1102
1103         ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold  ) ) {
1104
1105       
1106
1107       Int_t iDigitInCluster = 0 ; 
1108
1109       
1110
1111       if  ( IsInTower(digit) ) {   
1112
1113         // start a new Tower RecPoint
1114
1115         if(fNumberOfTowerClusters >= towerRecPoints->GetSize()) 
1116
1117           towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ;
1118
1119         
1120
1121         towerRecPoints->AddAt(new  AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ;
1122
1123         clu = dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(fNumberOfTowerClusters)) ; 
1124
1125         fNumberOfTowerClusters++ ; 
1126
1127         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ; 
1128
1129         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
1130
1131         iDigitInCluster++ ; 
1132
1133         digitsC->Remove(digit) ; 
1134
1135         
1136
1137       } else { 
1138
1139         
1140
1141         // start a new Pre Shower cluster
1142
1143         if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize()) 
1144
1145           preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1);
1146
1147         
1148
1149         preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ;
1150
1151         
1152
1153         clu =  dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(fNumberOfPreShoClusters))  ;  
1154
1155         fNumberOfPreShoClusters++ ; 
1156
1157         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) );     
1158
1159         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
1160
1161         iDigitInCluster++ ; 
1162
1163         digitsC->Remove(digit) ; 
1164
1165         nextdigit.Reset() ;
1166
1167         
1168
1169         // Here we remove remaining Tower digits, which cannot make a cluster
1170
1171         
1172
1173         if( notremoved ) { 
1174
1175           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
1176
1177             if( IsInTower(digit) )
1178
1179               digitsC->Remove(digit) ;
1180
1181             else 
1182
1183               break ; 
1184
1185           }
1186
1187           notremoved = kFALSE ;
1188
1189         }
1190
1191         
1192
1193       } // else        
1194
1195       
1196
1197       nextdigit.Reset() ;
1198
1199       
1200
1201       AliEMCALDigit * digitN ; 
1202
1203       index = 0 ;
1204
1205       while (index < iDigitInCluster){ // scan over digits already in cluster 
1206
1207         digit =  (AliEMCALDigit*)digits->At(clusterdigitslist[index])  ;      
1208
1209         index++ ; 
1210
1211         while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits 
1212
1213           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
1214
1215          switch (ineb ) {
1216
1217           case 0 :   // not a neighbour
1218
1219             break ;
1220
1221           case 1 :   // are neighbours 
1222
1223             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ;
1224
1225             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
1226
1227             iDigitInCluster++ ; 
1228
1229             digitsC->Remove(digitN) ;
1230
1231             break ;
1232
1233           case 2 :   // too far from each other
1234
1235             goto endofloop;   
1236
1237           } // switch
1238
1239           
1240
1241         } // while digitN
1242
1243         
1244
1245       endofloop: ;
1246
1247         nextdigit.Reset() ; 
1248
1249         
1250
1251       } // loop over cluster     
1252
1253       
1254
1255     } // energy theshold  
1256
1257     
1258
1259     
1260
1261   } // while digit
1262
1263   
1264
1265   delete digitsC ;
1266
1267   
1268
1269 }
1270
1271
1272
1273 //____________________________________________________________________________
1274
1275 void AliEMCALClusterizerv1::MakeUnfolding()
1276
1277 {
1278
1279   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
1280
1281   
1282
1283 //   // Unfolds clusters using the shape of an ElectroMagnetic shower
1284
1285 //   // Performs unfolding of all EMC/CPV clusters
1286
1287
1288
1289 //   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
1290
1291   
1292
1293 //   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
1294
1295 //   TObjArray * emcRecPoints = gime->TowerRecPoints() ; 
1296
1297 //   TObjArray * cpvRecPoints = gime->PreShoRecPoints() ; 
1298
1299 //   TClonesArray * digits = gime->Digits() ; 
1300
1301   
1302
1303 //   // Unfold first EMC clusters 
1304
1305 //   if(fNumberOfTowerClusters > 0){
1306
1307
1308
1309 //     Int_t nModulesToUnfold = geom->GetNModules() ; 
1310
1311
1312
1313 //     Int_t numberofNotUnfolded = fNumberOfTowerClusters ; 
1314
1315 //     Int_t index ;   
1316
1317 //     for(index = 0 ; index < numberofNotUnfolded ; index++){
1318
1319       
1320
1321 //       AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint *) emcRecPoints->At(index) ;
1322
1323 //       if(emcRecPoint->GetEMCALMod()> nModulesToUnfold)
1324
1325 //      break ;
1326
1327       
1328
1329 //       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
1330
1331 //       Int_t * maxAt = new Int_t[nMultipl] ;
1332
1333 //       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
1334
1335 //       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fTowerLocMaxCut,digits) ;
1336
1337       
1338
1339 //       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
1340
1341 //      UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
1342
1343 //      emcRecPoints->Remove(emcRecPoint); 
1344
1345 //      emcRecPoints->Compress() ;
1346
1347 //      index-- ;
1348
1349 //      fNumberOfTowerClusters -- ;
1350
1351 //      numberofNotUnfolded-- ;
1352
1353 //       }
1354
1355       
1356
1357 //       delete[] maxAt ; 
1358
1359 //       delete[] maxAtEnergy ; 
1360
1361 //     }
1362
1363 //   } 
1364
1365 //   // Unfolding of EMC clusters finished
1366
1367
1368
1369
1370
1371 //   // Unfold now CPV clusters
1372
1373 //   if(fNumberOfPreShoClusters > 0){
1374
1375     
1376
1377 //     Int_t nModulesToUnfold = geom->GetNModules() ;
1378
1379
1380
1381 //     Int_t numberofPreShoNotUnfolded = fNumberOfPreShoClusters ;     
1382
1383 //     Int_t index ;   
1384
1385 //     for(index = 0 ; index < numberofPreShoNotUnfolded ; index++){
1386
1387       
1388
1389 //       AliEMCALRecPoint * recPoint = (AliEMCALRecPoint *) cpvRecPoints->At(index) ;
1390
1391
1392
1393 //       if(recPoint->GetEMCALMod()> nModulesToUnfold)
1394
1395 //      break ;
1396
1397       
1398
1399 //       AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint*) recPoint ; 
1400
1401       
1402
1403 //       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
1404
1405 //       Int_t * maxAt = new Int_t[nMultipl] ;
1406
1407 //       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
1408
1409 //       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fPreShoLocMaxCut,digits) ;
1410
1411       
1412
1413 //       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
1414
1415 //      UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
1416
1417 //      cpvRecPoints->Remove(emcRecPoint); 
1418
1419 //      cpvRecPoints->Compress() ;
1420
1421 //      index-- ;
1422
1423 //      numberofPreShoNotUnfolded-- ;
1424
1425 //      fNumberOfPreShoClusters-- ;
1426
1427 //       }
1428
1429       
1430
1431 //       delete[] maxAt ; 
1432
1433 //       delete[] maxAtEnergy ; 
1434
1435 //     } 
1436
1437 //   }
1438
1439 //   //Unfolding of PreSho clusters finished
1440
1441   
1442
1443 }
1444
1445
1446
1447 //____________________________________________________________________________
1448
1449 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
1450
1451
1452
1453   // Shape of the shower (see EMCAL TDR)
1454
1455   // If you change this function, change also the gradient evaluation in ChiSquare()
1456
1457
1458
1459   Double_t r4    = r*r*r*r ;
1460
1461   Double_t r295  = TMath::Power(r, 2.95) ;
1462
1463   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
1464
1465   return shape ;
1466
1467 }
1468
1469
1470
1471 //____________________________________________________________________________
1472
1473 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, 
1474
1475                                                  Int_t nMax, 
1476
1477                                                  AliEMCALDigit ** maxAt, 
1478
1479                                                  Float_t * maxAtEnergy)
1480
1481 {
1482
1483   // Performs the unfolding of a cluster with nMax overlapping showers 
1484
1485   
1486
1487   Fatal("AliEMCALClusterizerv1::UnfoldCluster", "--> Unfolding not implemented") ;
1488
1489
1490
1491  //  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
1492
1493 //   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
1494
1495 //   const TClonesArray * digits = gime->Digits() ; 
1496
1497 //   TObjArray * emcRecPoints = gime->TowerRecPoints() ; 
1498
1499 //   TObjArray * cpvRecPoints = gime->PreShoRecPoints() ; 
1500
1501
1502
1503 //   Int_t nPar = 3 * nMax ;
1504
1505 //   Float_t * fitparameters = new Float_t[nPar] ;
1506
1507
1508
1509 //   Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
1510
1511 //   if( !rv ) {
1512
1513 //     // Fit failed, return and remove cluster
1514
1515 //     delete[] fitparameters ; 
1516
1517 //     return ;
1518
1519 //   }
1520
1521
1522
1523 //   // create ufolded rec points and fill them with new energy lists
1524
1525 //   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
1526
1527 //   // and later correct this number in acordance with actual energy deposition
1528
1529
1530
1531 //   Int_t nDigits = iniTower->GetMultiplicity() ;  
1532
1533 //   Float_t * efit = new Float_t[nDigits] ;
1534
1535 //   Float_t xDigit=0.,zDigit=0.,distance=0. ;
1536
1537 //   Float_t xpar=0.,zpar=0.,epar=0.  ;
1538
1539 //   Int_t relid[4] ;
1540
1541 //   AliEMCALDigit * digit = 0 ;
1542
1543 //   Int_t * emcDigits = iniTower->GetDigitsList() ;
1544
1545
1546
1547 //   Int_t iparam ;
1548
1549 //   Int_t iDigit ;
1550
1551 //   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1552
1553 //     digit = (AliEMCALDigit*) digits->At(emcDigits[iDigit] ) ;   
1554
1555 //     geom->AbsToRelNumbering(digit->GetId(), relid) ;
1556
1557 //     geom->RelPosInModule(relid, xDigit, zDigit) ;
1558
1559 //     efit[iDigit] = 0;
1560
1561
1562
1563 //     iparam = 0 ;    
1564
1565 //     while(iparam < nPar ){
1566
1567 //       xpar = fitparameters[iparam] ;
1568
1569 //       zpar = fitparameters[iparam+1] ;
1570
1571 //       epar = fitparameters[iparam+2] ;
1572
1573 //       iparam += 3 ;
1574
1575 //       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1576
1577 //       distance =  TMath::Sqrt(distance) ;
1578
1579 //       efit[iDigit] += epar * ShowerShape(distance) ;
1580
1581 //     }
1582
1583 //   }
1584
1585   
1586
1587
1588
1589 //   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
1590
1591 //   // so that energy deposited in each cell is distributed betwin new clusters proportionally
1592
1593 //   // to its contribution to efit
1594
1595
1596
1597 //   Float_t * emcEnergies = iniTower->GetEnergiesList() ;
1598
1599 //   Float_t ratio ;
1600
1601
1602
1603 //   iparam = 0 ;
1604
1605 //   while(iparam < nPar ){
1606
1607 //     xpar = fitparameters[iparam] ;
1608
1609 //     zpar = fitparameters[iparam+1] ;
1610
1611 //     epar = fitparameters[iparam+2] ;
1612
1613 //     iparam += 3 ;    
1614
1615     
1616
1617 //     AliEMCALTowerRecPoint * emcRP = 0 ;  
1618
1619
1620
1621 //     if(iniTower->IsTower()){ //create new entries in fTowerRecPoints...
1622
1623       
1624
1625 //       if(fNumberOfTowerClusters >= emcRecPoints->GetSize())
1626
1627 //      emcRecPoints->Expand(2*fNumberOfTowerClusters) ;
1628
1629       
1630
1631 //       (*emcRecPoints)[fNumberOfTowerClusters] = new AliEMCALTowerRecPoint("") ;
1632
1633 //       emcRP = (AliEMCALTowerRecPoint *) emcRecPoints->At(fNumberOfTowerClusters);
1634
1635 //       fNumberOfTowerClusters++ ;
1636
1637 //     }
1638
1639 //     else{//create new entries in fPreShoRecPoints
1640
1641 //       if(fNumberOfPreShoClusters >= cpvRecPoints->GetSize())
1642
1643 //      cpvRecPoints->Expand(2*fNumberOfPreShoClusters) ;
1644
1645       
1646
1647 //       (*cpvRecPoints)[fNumberOfPreShoClusters] = new AliEMCALPreShoRecPoint("") ;
1648
1649 //       emcRP = (AliEMCALTowerRecPoint *) cpvRecPoints->At(fNumberOfPreShoClusters);
1650
1651 //       fNumberOfPreShoClusters++ ;
1652
1653 //     }
1654
1655     
1656
1657 //     Float_t eDigit ;
1658
1659 //     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1660
1661 //       digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ; 
1662
1663 //       geom->AbsToRelNumbering(digit->GetId(), relid) ;
1664
1665 //       geom->RelPosInModule(relid, xDigit, zDigit) ;
1666
1667 //       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1668
1669 //       distance =  TMath::Sqrt(distance) ;
1670
1671 //       ratio = epar * ShowerShape(distance) / efit[iDigit] ; 
1672
1673 //       eDigit = emcEnergies[iDigit] * ratio ;
1674
1675 //       emcRP->AddDigit( *digit, eDigit ) ;
1676
1677 //     }        
1678
1679 //   }
1680
1681  
1682
1683 //   delete[] fitparameters ; 
1684
1685 //   delete[] efit ; 
1686
1687   
1688
1689 }
1690
1691
1692
1693 //_____________________________________________________________________________
1694
1695 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1696
1697 {
1698
1699   // Calculates the Chi square for the cluster unfolding minimization
1700
1701   // Number of parameters, Gradient, Chi squared, parameters, what to do
1702
1703
1704
1705   abort() ; 
1706
1707  //  Fatal("AliEMCALClusterizerv1::UnfoldingChiSquare","-->Unfolding not implemented") ;
1708
1709
1710
1711 //   TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
1712
1713
1714
1715 //   AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint*) toMinuit->At(0)  ;
1716
1717 //   TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
1718
1719
1720
1721
1722
1723   
1724
1725 //   //  AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint *) gMinuit->GetObjectFit() ; // TowerRecPoint to fit
1726
1727
1728
1729 //   Int_t * emcDigits     = emcRP->GetDigitsList() ;
1730
1731
1732
1733 //   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ; 
1734
1735
1736
1737 //   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1738
1739
1740
1741 //   const AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ; 
1742
1743 //   fret = 0. ;     
1744
1745 //   Int_t iparam ;
1746
1747
1748
1749 //   if(iflag == 2)
1750
1751 //     for(iparam = 0 ; iparam < nPar ; iparam++)    
1752
1753 //       Grad[iparam] = 0 ; // Will evaluate gradient
1754
1755   
1756
1757 //   Double_t efit ;    
1758
1759
1760
1761 //   AliEMCALDigit * digit ;
1762
1763 //   Int_t iDigit ;
1764
1765
1766
1767 //   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1768
1769
1770
1771 //     digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ; 
1772
1773
1774
1775 //     Int_t relid[4] ;
1776
1777 //     Float_t xDigit ;
1778
1779 //     Float_t zDigit ;
1780
1781
1782
1783 //     geom->AbsToRelNumbering(digit->GetId(), relid) ;
1784
1785
1786
1787 //     geom->RelPosInModule(relid, xDigit, zDigit) ;
1788
1789
1790
1791 //      if(iflag == 2){  // calculate gradient
1792
1793 //        Int_t iParam = 0 ;
1794
1795 //        efit = 0 ;
1796
1797 //        while(iParam < nPar ){
1798
1799 //       Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
1800
1801 //       iParam++ ; 
1802
1803 //       distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
1804
1805 //       distance = TMath::Sqrt( distance ) ; 
1806
1807 //       iParam++ ;      
1808
1809 //       efit += x[iParam] * ShowerShape(distance) ;
1810
1811 //       iParam++ ;
1812
1813 //        }
1814
1815 //        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) 
1816
1817 //        iParam = 0 ;
1818
1819 //        while(iParam < nPar ){
1820
1821 //       Double_t xpar = x[iParam] ;
1822
1823 //       Double_t zpar = x[iParam+1] ;
1824
1825 //       Double_t epar = x[iParam+2] ;
1826
1827 //       Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1828
1829 //       Double_t shape = sum * ShowerShape(dr) ;
1830
1831 //       Double_t r4 = dr*dr*dr*dr ;
1832
1833 //       Double_t r295 = TMath::Power(dr,2.95) ;
1834
1835 //       Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1836
1837 //                                       0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1838
1839          
1840
1841 //       Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
1842
1843 //       iParam++ ; 
1844
1845 //       Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
1846
1847 //       iParam++ ; 
1848
1849 //       Grad[iParam] += shape ;                                  // Derivative over energy             
1850
1851 //       iParam++ ; 
1852
1853 //        }
1854
1855 //      }
1856
1857 //      efit = 0;
1858
1859 //      iparam = 0 ;
1860
1861
1862
1863 //      while(iparam < nPar ){
1864
1865 //        Double_t xpar = x[iparam] ;
1866
1867 //        Double_t zpar = x[iparam+1] ;
1868
1869 //        Double_t epar = x[iparam+2] ;
1870
1871 //        iparam += 3 ;
1872
1873 //        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1874
1875 //        distance =  TMath::Sqrt(distance) ;
1876
1877 //        efit += epar * ShowerShape(distance) ;
1878
1879 //      }
1880
1881
1882
1883 //      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
1884
1885 //      // Here we assume, that sigma = sqrt(E)
1886
1887 //   }
1888
1889
1890
1891 }
1892
1893
1894
1895 //____________________________________________________________________________
1896
1897 void AliEMCALClusterizerv1::Print(Option_t * option)const
1898
1899 {
1900
1901   // Print clusterizer parameters
1902
1903
1904
1905   if( strcmp(GetName(), "") !=0 ){
1906
1907     
1908
1909     // Print parameters
1910
1911  
1912
1913     TString taskName(GetName()) ; 
1914
1915     taskName.ReplaceAll(Version(), "") ;
1916
1917
1918
1919     cout << "---------------"<< taskName.Data() << " " << GetTitle()<< "-----------" << endl 
1920
1921          << "Clusterizing digits from the file: " << taskName.Data() << endl 
1922
1923          << "                           Branch: " << GetName() << endl 
1924
1925          << endl 
1926
1927          << "                       EMC Clustering threshold = " << fTowerClusteringThreshold << endl
1928
1929          << "                       EMC Local Maximum cut    = " << fTowerLocMaxCut << endl
1930
1931          << "                       EMC Logarothmic weight   = " << fW0 << endl
1932
1933          << endl
1934
1935          << "                       CPV Clustering threshold = " << fPreShoClusteringThreshold << endl
1936
1937          << "                       CPV Local Maximum cut    = " << fPreShoLocMaxCut << endl
1938
1939        << "                       CPV Logarothmic weight   = " << fW0CPV << endl
1940
1941          << endl ;
1942
1943     if(fToUnfold)
1944
1945       cout << " Unfolding on " << endl ;
1946
1947     else
1948
1949       cout << " Unfolding off " << endl ;
1950
1951     
1952
1953     cout << "------------------------------------------------------------------" <<endl ;
1954
1955   }
1956
1957   else
1958
1959     cout << " AliEMCALClusterizerv1 not initialized " << endl ;
1960
1961 }
1962
1963 //____________________________________________________________________________
1964
1965 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
1966
1967 {
1968
1969   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
1970
1971
1972
1973   TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ; 
1974
1975   TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ; 
1976
1977
1978
1979   cout << "AliEMCALClusterizerv1: : event "<<gAlice->GetEvNumber() << endl ;
1980
1981   cout << "       Found "<< towerRecPoints->GetEntriesFast() << " TOWER Rec Points and " 
1982
1983            << preshoRecPoints->GetEntriesFast() << " PRE SHOWER RecPoints" << endl ;
1984
1985
1986
1987   fRecPointsInRun +=  towerRecPoints->GetEntriesFast() ; 
1988
1989   fRecPointsInRun +=  preshoRecPoints->GetEntriesFast() ; 
1990
1991
1992
1993   if(strstr(option,"all")) {
1994
1995
1996
1997     cout << "Tower clusters " << endl ;
1998
1999     cout << " Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list "      <<  endl;      
2000
2001     
2002
2003     Int_t index ;
2004
2005     for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) {
2006
2007       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(towerRecPoints->At(index)) ; 
2008
2009       TVector3  globalpos;  
2010
2011       rp->GetGlobalPosition(globalpos);
2012
2013       Float_t lambda[2]; 
2014
2015       rp->GetElipsAxis(lambda);
2016
2017       Int_t * primaries; 
2018
2019       Int_t nprimaries;
2020
2021       primaries = rp->GetPrimaries(nprimaries);
2022
2023
2024
2025       cout << setw(4) << rp->GetIndexInList() << "   " 
2026
2027            << setw(7) << setprecision(3) << rp->GetEnergy() << "      "
2028
2029            << setw(3) << rp->GetMultiplicity() << "      " 
2030
2031            << setw(1) << rp->GetEMCALArm() << "     " 
2032
2033            << setw(5) << setprecision(4) << globalpos.X() << "  " 
2034
2035            << setw(5) << setprecision(4) << globalpos.Y() << "  " 
2036
2037            << setw(5) << setprecision(4) << globalpos.Z() << "     "
2038
2039            << setw(4) << setprecision(2) << lambda[0]  << "  "
2040
2041            << setw(4) << setprecision(2) << lambda[1]  << "  "
2042
2043            << setw(2) << nprimaries << "  " ;
2044
2045      
2046
2047       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
2048
2049         cout << setw(4) <<   primaries[iprimary] << "  "  ;
2050
2051       cout << endl ;     
2052
2053     }
2054
2055
2056
2057     //Now plot Pre shower recPoints
2058
2059
2060
2061     cout << "-----------------------------------------------------------------------"<<endl ;
2062
2063
2064
2065     cout << "PreShower clusters " << endl ;
2066
2067     cout << " Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list "      <<  endl;      
2068
2069     
2070
2071     for (index = 0 ; index < preshoRecPoints->GetEntries() ; index++) {
2072
2073       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(index)) ; 
2074
2075       TVector3  globalpos;  
2076
2077       rp->GetGlobalPosition(globalpos);
2078
2079       Float_t lambda[2]; 
2080
2081       rp->GetElipsAxis(lambda);
2082
2083       Int_t * primaries; 
2084
2085       Int_t nprimaries;
2086
2087       primaries = rp->GetPrimaries(nprimaries);
2088
2089
2090
2091       cout << setw(4) << rp->GetIndexInList() << "   " 
2092
2093            << setw(7) << setprecision(3) << rp->GetEnergy() << "      "
2094
2095            << setw(3) << rp->GetMultiplicity() << "      " 
2096
2097            << setw(1) << rp->GetEMCALArm() << "     " 
2098
2099            << setw(5) << setprecision(4) << globalpos.X() << "  " 
2100
2101            << setw(5) << setprecision(4) << globalpos.Y() << "  " 
2102
2103            << setw(5) << setprecision(4) << globalpos.Z() << "     "
2104
2105            << setw(4) << setprecision(2) << lambda[0]  << "  "
2106
2107            << setw(4) << setprecision(2) << lambda[1]  << "  "
2108
2109            << setw(2) << nprimaries << "  " ;
2110
2111      
2112
2113       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
2114
2115         cout << setw(4) <<   primaries[iprimary] << "  "  ;
2116
2117       cout << endl ;     
2118
2119     }
2120
2121
2122
2123     cout << "-----------------------------------------------------------------------"<<endl ;
2124
2125   }
2126
2127 }
2128
2129
2130