2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
19 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
20 // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
21 // of new IO (à la PHOS)
22 //////////////////////////////////////////////////////////////////////////////
23 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
24 // unfolds the clusters having several local maxima.
25 // Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
26 // EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
27 // parameters including input digits branch title, thresholds etc.)
28 // This TTask is normally called from Reconstructioner, but can as well be used in
31 // root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")
32 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
33 // //reads gAlice from header file "..."
34 // root [1] cl->ExecuteTask()
35 // //finds RecPoints in all events stored in galice.root
36 // root [2] cl->SetDigitsBranch("digits2")
37 // //sets another title for Digitis (input) branch
38 // root [3] cl->SetRecPointsBranch("recp2")
39 // //sets another title four output branches
40 // root [4] cl->SetTowerLocalMaxCut(0.03)
41 // //set clusterization parameters
42 // root [5] cl->ExecuteTask("deb all time")
43 // //once more finds RecPoints options are
44 // // deb - print number of found rec points
45 // // deb all - print number of found RecPoints and some their characteristics
46 // // time - print benchmarking results
48 // --- ROOT system ---
57 #include "TBenchmark.h"
59 // --- Standard library ---
62 // --- AliRoot header files ---
63 #include "AliEMCALGetter.h"
64 #include "AliEMCALClusterizerv1.h"
65 #include "AliEMCALRecPoint.h"
66 #include "AliEMCALDigit.h"
67 #include "AliEMCALDigitizer.h"
69 #include "AliEMCALGeometry.h"
71 ClassImp(AliEMCALClusterizerv1)
73 Int_t addOn[20][60][60];
75 //____________________________________________________________________________
76 AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
78 // default ctor (to be used mainly by Streamer)
81 fDefaultInit = kTRUE ;
82 for(Int_t is=0;is<20;is++){
83 for(Int_t i=0;i<60;i++){
84 for(Int_t j=0;j<60;j++){
89 //PH cout<<"file to read 1"<<endl;
91 //PH cout<<"file read 1"<<endl;
94 //____________________________________________________________________________
95 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
96 :AliEMCALClusterizer(alirunFileName, eventFolderName)
98 // ctor with the indication of the file where header Tree and digits Tree are stored
102 fDefaultInit = kFALSE ;
103 for(Int_t is=0;is<20;is++){
104 for(Int_t i=0;i<60;i++){
105 for(Int_t j=0;j<60;j++){
110 //PH cout<<"file to read 2"<<endl;
112 //PH cout<<"file read 2"<<endl;
116 //____________________________________________________________________________
117 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
123 //____________________________________________________________________________
124 const TString AliEMCALClusterizerv1::BranchName() const
130 //____________________________________________________________________________
131 Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp) const
133 //To be replased later by the method, reading individual parameters from the database
134 return -fADCpedestalECA + amp * fADCchannelECA ;
137 //____________________________________________________________________________
138 void AliEMCALClusterizerv1::Exec(Option_t * option)
140 // Steering method to perform clusterization for events
141 // in the range from fFirstEvent to fLastEvent.
142 // This range is optionally set by SetEventRange().
143 // if fLastEvent=-1 (by default), then process events until the end.
145 if(strstr(option,"tim"))
146 gBenchmark->Start("EMCALClusterizer");
148 if(strstr(option,"print"))
151 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
153 if (fLastEvent == -1)
154 fLastEvent = gime->MaxEvent() - 1;
155 Int_t nEvents = fLastEvent - fFirstEvent + 1;
158 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
159 gime->Event(ievent,"D") ;
160 GetCalibrationParameters() ;
162 fNumberOfECAClusters = 0;
171 if(strstr(option,"deb"))
172 PrintRecPoints(option) ;
174 //increment the total number of recpoints per run
175 fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ;
180 if(strstr(option,"tim")){
181 gBenchmark->Stop("EMCALClusterizer");
182 printf("Exec took %f seconds for Clusterizing %f seconds per event",
183 gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ;
190 //____________________________________________________________________________
191 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
192 Int_t nPar, Float_t * fitparameters) const
194 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
195 // The initial values for fitting procedure are set equal to the positions of local maxima.
196 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
198 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
199 TClonesArray * digits = gime->Digits() ;
201 gMinuit->mncler(); // Reset Minuit's list of paramters
202 gMinuit->SetPrintLevel(-1) ; // No Printout
203 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
204 // To set the address of the minimization function
205 TList * toMinuit = new TList();
206 toMinuit->AddAt(emcRP,0) ;
207 toMinuit->AddAt(digits,1) ;
209 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
211 // filling initial values for fit parameters
212 AliEMCALDigit * digit ;
216 Int_t nDigits = (Int_t) nPar / 3 ;
220 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
222 for(iDigit = 0; iDigit < nDigits; iDigit++){
223 digit = maxAt[iDigit];
228 geom->AbsToRelNumbering(digit->GetId(), relid) ;
229 geom->PosInAlice(relid, x, z) ;
231 Float_t energy = maxAtEnergy[iDigit] ;
233 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
236 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ;
239 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
242 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ;
245 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
248 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;
253 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
258 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
259 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
260 gMinuit->SetMaxIterations(5);
261 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
263 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
265 if(ierflg == 4){ // Minimum not found
266 Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
269 for(index = 0; index < nPar; index++){
272 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
273 fitparameters[index] = val ;
281 //____________________________________________________________________________
282 void AliEMCALClusterizerv1::GetCalibrationParameters()
284 // Gets the parameters for the calibration from the digitizer
285 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
287 if ( !gime->Digitizer() )
288 gime->LoadDigitizer();
289 AliEMCALDigitizer * dig = gime->Digitizer();
291 fADCchannelECA = dig->GetECAchannel() ;
292 fADCpedestalECA = dig->GetECApedestal();
293 //PH cout<<"ChannelECA, peds "<<fADCchannelECA<<" "<<fADCpedestalECA<<endl;
296 //____________________________________________________________________________
297 void AliEMCALClusterizerv1::Init()
299 // Make all memory allocations which can not be done in default constructor.
300 // Attach the Clusterizer task to the list of EMCAL tasks
302 AliEMCALGetter * gime = AliEMCALGetter::Instance();
304 gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
306 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
307 //PH cout<<"gime,geom "<<gime<<" "<<geom<<endl;
309 //Sub fNTowers = geom->GetNZ() * geom->GetNPhi() ;
312 gMinuit = new TMinuit(100) ;
313 //Sub if ( !gime->Clusterizer() )
314 //Sub gime->PostClusterizer(this);
316 //PH cout<<"hists booked "<<endl;
319 //____________________________________________________________________________
320 void AliEMCALClusterizerv1::InitParameters()
322 // Initializes the parameters for the Clusterizer
323 fNumberOfECAClusters = 0;
325 fECAClusteringThreshold = 0.5; // value obtained from Alexei
326 fECALocMaxCut = 0.03;
331 fRecPointsInRun = 0 ;
335 //____________________________________________________________________________
336 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
338 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
340 // = 2 are not neighbour but do not continue searching
341 // neighbours are defined as digits having at least a common vertex
342 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
343 // which is compared to a digit (d2) not yet in a cluster
345 AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
350 //Sub geom->AbsToRelNumbering(d1->GetId(), relid1) ;
357 geom->GetCellIndex(d1->GetId(), nSupMod,nTower,nIphi,nIeta);
358 geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
369 geom->GetCellIndex(d2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
370 geom->GetCellPhiEtaIndexInSModule(nSupMod1, nTower1,nIphi1,nIeta1, iphi1,ieta1);
375 //Sub geom->AbsToRelNumbering(d2->GetId(), relid2) ;
378 Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
379 Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
381 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
385 if((relid2[0] > relid1[0]) && (relid2[1] > relid1[1]+1))
386 rv = 2; // Difference in row numbers is too large to look further
390 if(rv==1)printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
391 rv, d1->GetId(), relid1[0], relid1[1],
392 d2->GetId(), relid2[0], relid2[1]) ;
397 //____________________________________________________________________________
398 void AliEMCALClusterizerv1::Unload()
400 // Unloads the Digits and RecPoints
401 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
402 gime->EmcalLoader()->UnloadDigits() ;
403 gime->EmcalLoader()->UnloadRecPoints() ;
406 //____________________________________________________________________________
407 void AliEMCALClusterizerv1::WriteRecPoints()
410 // Creates new branches with given title
411 // fills and writes into TreeR.
413 AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
415 TObjArray * aECARecPoints = gime->ECARecPoints() ;
417 TClonesArray * digits = gime->Digits() ;
418 TTree * treeR = gime->TreeR(); ;
422 //Evaluate position, dispersion and other RecPoint properties for EC section
423 for(index = 0; index < aECARecPoints->GetEntries(); index++)
424 (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
426 aECARecPoints->Sort() ;
428 for(index = 0; index < aECARecPoints->GetEntries(); index++)
429 (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
431 aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ;
433 Int_t bufferSize = 32000 ;
434 Int_t splitlevel = 0 ;
437 TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
438 branchECA->SetTitle(BranchName());
442 gime->WriteRecPoints("OVERWRITE");
443 gime->WriteClusterizer("OVERWRITE");
446 //____________________________________________________________________________
447 void AliEMCALClusterizerv1::MakeClusters()
449 // Steering method to construct the clusters stored in a list of Reconstructed Points
450 // A cluster is defined as a list of neighbour digits
452 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
454 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
456 TObjArray * aECARecPoints = gime->ECARecPoints() ;
458 aECARecPoints->Delete() ;
460 TClonesArray * digits = gime->Digits() ;
461 TClonesArray * digitsC = dynamic_cast<TClonesArray*>(digits->Clone()) ;
463 // Clusterization starts
464 TIter nextdigit(digitsC) ;
465 AliEMCALDigit * digit;
469 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
470 digitAmp->Fill(digit->GetAmp());
475 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
476 AliEMCALRecPoint * clu = 0 ;
478 TArrayI clusterECAdigitslist(50000);
480 //Sub if (gDebug == 2) {
481 // printf("MakeClusters: id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()),
482 // fECAClusteringThreshold) ;
484 ////////////////////////temp solution, embedding///////////////////
485 int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
487 geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
488 geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
490 ///////////////////////////////////
492 //cout<<ieta<<" "<<iphi<<endl;
494 // if ( geom->IsInECA(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold ) ){
495 if (geom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) > fECAClusteringThreshold ) ){
496 //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"11 digit, add "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
497 // cout<<"crossed the threshold "<<endl;
498 Int_t iDigitInECACluster = 0;
499 // start a new Tower RecPoint
500 if(fNumberOfECAClusters >= aECARecPoints->GetSize())
501 aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
502 AliEMCALRecPoint * rp = new AliEMCALRecPoint("") ;
503 aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
504 clu = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
505 fNumberOfECAClusters++ ;
506 clu->AddDigit(*digit, Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1])) ;
507 clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
508 iDigitInECACluster++ ;
509 // cout<<" 1st setp:cluno, digNo "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
510 digitsC->Remove(digit) ;
512 printf("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ;
515 AliEMCALDigit * digitN ;
518 // Find the neighbours
519 while (index < iDigitInECACluster){ // scan over digits already in cluster
520 digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]) ;
522 while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits
523 // cout<<"we have new digit "<<endl;
524 // check that the digit is above the min E Cut
526 geom->GetCellIndex(digitN->GetId(), nSupMod,nTower,nIphi,nIeta);
527 geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
529 if(Calibrate(digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut ) digitsC->Remove(digitN);
530 // cout<<" new digit above ECut "<<endl;
531 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
532 // cout<<" new digit neighbour?? "<<ineb<<endl;
534 case 0 : // not a neighbour
536 case 1 : // are neighbours
537 //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"22 digit, add "<<nSupMod<<" "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
538 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) ) ;
539 clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
540 iDigitInECACluster++ ;
541 // cout<<"when neighbour: cluno, digNo "<<digit->GetId()<<" "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
542 digitsC->Remove(digitN) ;
544 // case 2 : // too far from each other
545 //Subh goto endofloop1;
546 // cout<<"earlier go to end of loop"<<endl;
548 //cout<<"in nextDigit loop "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
553 } // loop over ECA cluster
554 } // energy threshold
555 else if(Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut ){
556 //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"33 digit, add "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
557 digitsC->Remove(digit);
559 //cout<<"after endofloop: cluno, digNo "<<fNumberOfECAClusters<<endl;
562 cout<<"total no of clusters "<<fNumberOfECAClusters<<" from "<<digits->GetEntriesFast()<<" digits"<<endl;
565 //____________________________________________________________________________
566 void AliEMCALClusterizerv1::MakeUnfolding() const
568 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
572 //____________________________________________________________________________
573 Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
575 // Shape of the shower (see EMCAL TDR)
576 // If you change this function, change also the gradient evaluation in ChiSquare()
578 Double_t r4 = r*r*r*r ;
579 Double_t r295 = TMath::Power(r, 2.95) ;
580 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
584 //____________________________________________________________________________
585 void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
587 AliEMCALDigit ** /*maxAt*/,
588 Float_t * /*maxAtEnergy*/) const
590 // Performs the unfolding of a cluster with nMax overlapping showers
592 Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
596 //_____________________________________________________________________________
597 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
599 Double_t * /*x*/, Int_t /*iflag*/)
601 // Calculates the Chi square for the cluster unfolding minimization
602 // Number of parameters, Gradient, Chi squared, parameters, what to do
604 ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
606 //____________________________________________________________________________
607 void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
609 // Print clusterizer parameters
611 TString message("\n") ;
613 if( strcmp(GetName(), "") !=0 ){
617 TString taskName(GetName()) ;
618 taskName.ReplaceAll(Version(), "") ;
620 printf("--------------- ");
621 printf(taskName.Data()) ;
624 printf("-----------\n");
625 printf("Clusterizing digits from the file: ");
626 printf(taskName.Data());
627 printf("\n Branch: ");
629 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
630 printf("\n ECA Logarithmic weight = %f", fECAW0);
632 printf("\nUnfolding on\n");
634 printf("\nUnfolding off\n");
636 printf("------------------------------------------------------------------");
639 printf("AliEMCALClusterizerv1 not initialized ") ;
642 //____________________________________________________________________________
643 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
645 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
647 TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ;
648 printf("PrintRecPoints: Clusterization result:") ;
650 printf("event # %d\n", gAlice->GetEvNumber() ) ;
651 printf(" Found %d ECA Rec Points\n ",
652 aECARecPoints->GetEntriesFast()) ;
654 fRecPointsInRun += aECARecPoints->GetEntriesFast() ;
656 if(strstr(option,"all")) {
658 printf("\n-----------------------------------------------------------------------\n") ;
659 printf("Clusters in ECAL section\n") ;
660 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
666 for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
667 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ;
669 //rp->GetGlobalPosition(globalpos);
671 rp->GetLocalPosition(localpos);
673 rp->GetElipsAxis(lambda);
676 primaries = rp->GetPrimaries(nprimaries);
677 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
678 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
679 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
680 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
682 if(rp->GetEnergy()>maxE){
683 maxE=rp->GetEnergy();
686 maxDis=rp->GetDispersion();
688 pointE->Fill(rp->GetEnergy());
689 pointL1->Fill(lambda[0]);
690 pointL2->Fill(lambda[1]);
691 pointDis->Fill(rp->GetDispersion());
692 pointMult->Fill(rp->GetMultiplicity());
694 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
695 printf("%d ", primaries[iprimary] ) ;
702 MaxDis->Fill(maxDis);
705 printf("\n-----------------------------------------------------------------------\n");
708 void AliEMCALClusterizerv1::BookHists()
710 pointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
711 pointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
712 pointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
713 pointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
714 pointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
715 digitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
716 MaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
717 MaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
718 MaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
719 MaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.);
721 void AliEMCALClusterizerv1::SaveHists()
723 recofile=new TFile("reco.root","RECREATE");
737 void AliEMCALClusterizerv1::ReadFile()
740 FILE *fp = fopen("hijing1.dat","r");
741 for(Int_t line=0;line<9113;line++){
744 ncols0 = fscanf(fp,"%d %d %d %d",&sm,&l1,&l2,&eg);
745 // cout<<eg<<" "<<l1<<" "<<l2<<endl;
746 addOn[sm-1][l1-1][l2-1]=eg;
747 //addOn[sm-1][l1-1][l2-1]=0;