[u/mrichter/AliRoot.git] / EMCAL / macros / TestESD.C


#if !defined(__CINT__) || defined(__MAKECINT__)

//Root include files 
//#include <Riostream.h>
#include <TFile.h>
//#include <TSystem.h>
#include <TH1F.h>
#include <TParticle.h>
#include <TRefArray.h>

//AliRoot include files 
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliESDCaloTrigger.h"
#include "AliPID.h"
#include "AliEMCALGeometry.h"

#endif

//Change the bool depending on what information you want to print
// when all FALSE, prints minimum cluster information.
Bool_t kPrintKine         = kFALSE; //Do not use for raw data.
Bool_t kPrintCaloCells    = kTRUE;
Bool_t kPrintCaloTrigger  = kFALSE;
Bool_t kPrintTrackMatches = kFALSE;
Bool_t kPrintClusterCells = kFALSE;
Bool_t kPrintClusterPID   = kFALSE;

void TestESD() {
  // Main method to read information stored in AliESDCaloClusters and AliESDCaloCells
	
  // Init some example histograms
  // ESD
  TH1F * hEta  = (TH1F*) new TH1F("hEta","reco eta",1000, -0.71,0.71); 
  TH1F * hPhi  = (TH1F*) new TH1F("hPhi","reco phi",130, 70,200); 
  TH1F * hE    = (TH1F*) new TH1F("hE"  ,"reco e",300, 0,30); 
  TH1F * hTime = (TH1F*) new TH1F("hTime"  ,"reco time",1000, 0,1000); 
  hEta ->SetXTitle("#eta");
  hPhi ->SetXTitle("#phi (deg)");
  hE   ->SetXTitle("E (GeV)");
  hTime->SetXTitle("time (ns)");
  
  // Monte Carlo
  TH1F * hMCEta = (TH1F*) new TH1F("hMCEta","MC eta",1000, -0.71,0.71); 
  TH1F * hMCPhi = (TH1F*) new TH1F("hMCPhi","MC phi",130, 70,200); 
  TH1F * hMCE   = (TH1F*) new TH1F("hMCE"  ,"MC e",300, 0,30); 
  hMCEta->SetXTitle("#eta");
  hMCPhi->SetXTitle("#phi (deg)");
  hMCE  ->SetXTitle("E (GeV)");
  
  //ESD - MonteCarlo
  TH1F * hDEta = (TH1F*) new TH1F("hDEta"," eta cluster - eta MC",500, -0.05,0.05); 
  TH1F * hDPhi = (TH1F*) new TH1F("hDPhi","phi cluster - phi MC",200, -20,20); 
  TH1F * hDE   = (TH1F*) new TH1F("hDE"  ,"e cluster - e MC",200, -10,10); 
  hDEta->SetXTitle("#eta_{reco}-#eta_{MC}");
  hDPhi->SetXTitle("#phi_{reco}-#phi_{MC} (deg)");
  hDE  ->SetXTitle("E_{reco}-E_{MC} (GeV)");
  
  //Open the ESD file, get the tree with events
  TFile* f = new TFile("AliESDs.root");
  TTree* esdTree = (TTree*)f->Get("esdTree");
  AliESDEvent* esd = new AliESDEvent();
  esd->ReadFromTree(esdTree);
  
  //Init geometry and array that will contain the clusters
  TRefArray* caloClusters = new TRefArray();
  AliEMCALGeometry *geom =  AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEAR") ;  	
  Float_t pos[3];
  
  //Loop of events
  Int_t nEvt = esdTree->GetEntries();
  for(Int_t iev = 0; iev < nEvt; iev++) {
    cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
    esdTree->GetEvent(iev);
    
    //Pass the geometry transformation matrix from ESDs to geometry
    for(Int_t mod=0; mod < (geom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
      //printf("matrix %d\n",mod);
      if(esd->GetEMCALMatrix(mod)) {
        //printf("EMCAL: mod %d, matrix %p\n",mod, esd->GetEMCALMatrix(mod));
        //(esd->GetEMCALMatrix(mod))->Print("");
        geom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
      }//matrix
    }//module
    
    //In case you want to play with MC data, get stack
    AliStack *stack = 0;
    if(kPrintKine){  
      AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),  "read");
      rl->LoadKinematics();  
      rl->GetEvent(iev);
      stack = rl->Stack();
    }  
    
    //get reconstructed vertex position
    Double_t vertex_position[3];
    esd->GetVertex()->GetXYZ(vertex_position);
    
    //------------------------------------------------------
    //Get Cells Array, all cells in event, print cells info 
    //------------------------------------------------------ 
    AliVCaloCells &cells= *(esd->GetEMCALCells());
    
    if(kPrintCaloCells){  
      Int_t nTotalCells = cells.GetNumberOfCells() ;  
      //Int_t type        = cells.GetType();
      for (Int_t icell=  0; icell <  nTotalCells; icell++) {
        cout<<"Cell   : "<<icell<<"/"<<nTotalCells<<" - ID: "<<cells.GetCellNumber(icell)<<"; Amplitude: "<<cells.GetAmplitude(icell)<<"; Time: "<<cells.GetTime(icell)*1e9;
        cout << "; MC label "<<cells.GetMCLabel(icell)<<"; Embeded E fraction "<<cells.GetEFraction(icell);
        cout<<endl;	       }// cell loop
    }
    
    //------------------------------------------------------
    // Calo Trigger 
    //------------------------------------------------------
    
	  if(kPrintCaloTrigger)
	  {  
		  AliESDCaloTrigger& trg = *(esd->GetCaloTrigger("EMCAL"));
		  
		  trg.Reset();
		  while (trg.Next())
		  {
			  int posX, posY;
			  trg.GetPosition(posX, posY);
			  
			  if (posX > -1 && posY > -1) 
			  {
				  Int_t ts = 0;
				  trg.GetL1TimeSum(ts);
				  
				  cout << "Position: " << posX << " " << posY << " L1 amplitude: " << ts << endl;
			  }
		  }
	  }
    
    //------------------------------------------------------
    // Calo Clusters 
    //------------------------------------------------------
    
    //Get CaloClusters Array
    esd->GetEMCALClusters(caloClusters);
    
    //loop over clusters
    Int_t nclus = caloClusters->GetEntries();
    for (Int_t icl = 0; icl < nclus; icl++) {
      
      AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
      Float_t energy = clus->E();
      clus->GetPosition(pos);
      TVector3 vpos(pos[0],pos[1],pos[2]);
      
      //We can get a momentum TLorentzVector per cluster, corrected by the vertex position 
      //TLorentzVector p;
      //clus->GetMomentum(p,vertex_position);
      
      Double_t cphi = vpos.Phi();
      if(cphi < 0) cphi +=TMath::TwoPi();
      Double_t ceta = vpos.Eta();
      
      Int_t nMatched   = clus->GetNTracksMatched();
      Int_t trackIndex = clus->GetTrackMatchedIndex();
      Int_t nLabels    = clus->GetNLabels();
      Int_t labelIndex = clus->GetLabel();
      Int_t nCells     = clus->GetNCells();
      
      //Fill some histograms
      hEta->Fill(ceta);
      hPhi->Fill(cphi*TMath::RadToDeg());
      hE  ->Fill(energy);
      hTime->Fill(clus->GetTOF()*1e9);
      
      //Print basic cluster information
      cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << "; Phi: " 
      << cphi*TMath::RadToDeg() << "; Eta: " << ceta << "; NCells: " << nCells << " ;#Matches: " << nMatched 
      << "; Index: " << trackIndex << "; #Labels: " << nLabels << " Index: " 
      << labelIndex << "; Time "<<clus->GetTOF()*1e9<<" ns "<<endl;
      
      //Print primary info
      if(stack && kPrintKine) {
        if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
          TParticle * particle = stack->Particle(labelIndex);
          //Fill histograms with primary info
          hMCEta->Fill(particle->Eta());
          hMCPhi->Fill(particle->Phi()*TMath::RadToDeg());
          hMCE  ->Fill(particle->Energy());
          hDEta ->Fill(ceta-particle->Eta());
          hDPhi ->Fill(cphi*TMath::RadToDeg()-particle->Phi()*TMath::RadToDeg());
          hDE   ->Fill(energy-particle->Energy());
          //Print primary values
          cout<<"         More  contributing primary: "<<particle->GetName()<<"; with kinematics: "<<endl;
          cout<<" \t     Energy: "<<particle->Energy()<<"; Phi: "<<particle->Phi()*TMath::RadToDeg()<<"; Eta: "<<particle->Eta()<<endl;   
          for(Int_t i = 1; i < nLabels; i++){
            //particle = stack->Particle((((AliESDCaloCluster*)clus)->GetLabelsArray())->At(i));
            particle = stack->Particle((clus->GetLabels())[i]);
            //or Int_t *labels = clus->GetLabels();
            //particle = stack->Particle(labels[i]);
            cout<<"         Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
          }
        }
        else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
        else cout<<"Negative label!!!  : "<<labelIndex<<endl;
      } // play with stack
      
      // Matching results
      if(kPrintTrackMatches && trackIndex >= 0) {
        AliESDtrack* track = esd->GetTrack(trackIndex);
        Double_t tphi = track->GetOuterParam()->Phi();
        Double_t teta = track->GetOuterParam()->Eta();
        Double_t tmom = track->GetOuterParam()->P();
        cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
        
        Double_t deta = teta - ceta;
        Double_t dphi = tphi - cphi;
        if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
        if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
        Double_t dR = sqrt(dphi*dphi + deta*deta);
        
        Double_t pOverE = tmom/energy;
        
        if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
          cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
        }
      }// matching
      
      //Get PID weights and print them
      if(kPrintClusterPID){
        const Double_t *pid = clus->GetPID();
        printf("PID weights: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
               pid[AliVCluster::kPhoton],   pid[AliVCluster::kPi0],
               pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
               pid[AliVCluster::kPion],     pid[AliVCluster::kKaon],   pid[AliVCluster::kProton],
               pid[AliVCluster::kNeutron],  pid[AliVCluster::kKaon0]);
      }//PID
      
      //Get CaloCells of cluster and print their info, position.
      if(kPrintClusterCells){	
        UShort_t * index    = clus->GetCellsAbsId() ;
        Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
        Int_t sm = -1;
        for(Int_t i = 0; i < nCells ; i++){
          Int_t absId       =   index[i]; // or clus->GetCellNumber(i) ;
          Double_t ampFract =  fraction[i];
          Float_t amp       = cells.GetCellAmplitude(absId) ;
          Double_t time     = cells.GetCellTime(absId);
          cout<<"\t Cluster Cell: AbsID : "<< absId << " == "<<clus->GetCellAbsId(i) <<"; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time*1e9<<endl;
          //Geometry methods  
          Int_t iSupMod =  0 ;
          Int_t iTower  =  0 ;
          Int_t iIphi   =  0 ;
          Int_t iIeta   =  0 ;
          Int_t iphi    =  0 ;
          Int_t ieta    =  0 ;
          if(geom){
            geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
            //Gives SuperModule and Tower numbers
            geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
                                              iIphi, iIeta,iphi,ieta);
            //Gives label of cell in eta-phi position per each supermodule
            Float_t cellPhi = 0;
            Float_t cellEta = 0;
            geom->EtaPhiFromIndex(absId,cellEta,cellPhi);
            cout<< "                SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
            <<"; Phi "<<iIphi<<"; Index: Cell Eta "<<ieta<<"; Cell Phi "<<iphi
            <<"; Global: Cell Eta "<<cellEta<<"; Cell Phi "<<cellPhi*TMath::RadToDeg()<<endl;
            if(i==0) sm = iSupMod;
            else{
              if(sm!=iSupMod) printf("******CLUSTER SHARED BY 2 SuperModules!!!!\n");
            }	
          }// geometry on
        }// cluster cell loop
      }// print cell clusters
    } //cluster loop
  } // event loop
  
  
  //Write histograms in a file
  TFile * fhisto = (TFile*) new TFile("histos.root","recreate");
  hEta->Write();
  hPhi->Write();
  hE  ->Write();
  hTime->Write();
  if(kPrintKine){
    hMCEta->Write();
    hMCPhi->Write();
    hMCE  ->Write();
    hDEta->Write();
    hDPhi->Write();
    hDE  ->Write();
  }
  fhisto->Close();
}
Commit	Line	Data
c8fe2783	1
cd259404	2	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	3
	4	//Root include files
c8fe2783	5	//#include <Riostream.h>
cd259404	6	#include <TFile.h>
c8fe2783	7	//#include <TSystem.h>
cd259404	8	#include <TH1F.h>
cd259404	9	#include <TParticle.h>
cd259404	10	#include <TRefArray.h>
cd259404	11
	12	//AliRoot include files
	13	#include "AliRunLoader.h"
	14	#include "AliStack.h"
	15	#include "AliESDEvent.h"
	16	#include "AliESDVertex.h"
	17	#include "AliESDCaloCluster.h"
	18	#include "AliESDCaloCells.h"
754a26d5	19	#include "AliESDCaloTrigger.h"
cd259404	20	#include "AliPID.h"
c8fe2783	21	#include "AliEMCALGeometry.h"
cd259404	22
	23	#endif
	24
9480984d	25	//Change the bool depending on what information you want to print
9480984d	26	// when all FALSE, prints minimum cluster information.
c8fe2783	27	Bool_t kPrintKine = kFALSE; //Do not use for raw data.
77e93dc2	28	Bool_t kPrintCaloCells = kTRUE;
754a26d5	29	Bool_t kPrintCaloTrigger = kFALSE;
9480984d	30	Bool_t kPrintTrackMatches = kFALSE;
9480984d	31	Bool_t kPrintClusterCells = kFALSE;
c8fe2783	32	Bool_t kPrintClusterPID = kFALSE;
9480984d	33
cd259404	34	void TestESD() {
c8fe2783	35	// Main method to read information stored in AliESDCaloClusters and AliESDCaloCells
	36
	37	// Init some example histograms
	38	// ESD
	39	TH1F * hEta = (TH1F*) new TH1F("hEta","reco eta",1000, -0.71,0.71);
	40	TH1F * hPhi = (TH1F*) new TH1F("hPhi","reco phi",130, 70,200);
	41	TH1F * hE = (TH1F*) new TH1F("hE" ,"reco e",300, 0,30);
	42	TH1F * hTime = (TH1F*) new TH1F("hTime" ,"reco time",1000, 0,1000);
	43	hEta ->SetXTitle("#eta");
	44	hPhi ->SetXTitle("#phi (deg)");
	45	hE ->SetXTitle("E (GeV)");
	46	hTime->SetXTitle("time (ns)");
	47
	48	// Monte Carlo
	49	TH1F * hMCEta = (TH1F*) new TH1F("hMCEta","MC eta",1000, -0.71,0.71);
	50	TH1F * hMCPhi = (TH1F*) new TH1F("hMCPhi","MC phi",130, 70,200);
	51	TH1F * hMCE = (TH1F*) new TH1F("hMCE" ,"MC e",300, 0,30);
	52	hMCEta->SetXTitle("#eta");
	53	hMCPhi->SetXTitle("#phi (deg)");
	54	hMCE ->SetXTitle("E (GeV)");
	55
	56	//ESD - MonteCarlo
	57	TH1F * hDEta = (TH1F*) new TH1F("hDEta"," eta cluster - eta MC",500, -0.05,0.05);
	58	TH1F * hDPhi = (TH1F*) new TH1F("hDPhi","phi cluster - phi MC",200, -20,20);
	59	TH1F * hDE = (TH1F*) new TH1F("hDE" ,"e cluster - e MC",200, -10,10);
	60	hDEta->SetXTitle("#eta_{reco}-#eta_{MC}");
	61	hDPhi->SetXTitle("#phi_{reco}-#phi_{MC} (deg)");
	62	hDE ->SetXTitle("E_{reco}-E_{MC} (GeV)");
	63
	64	//Open the ESD file, get the tree with events
cd259404	65	TFile* f = new TFile("AliESDs.root");
cd259404	66	TTree* esdTree = (TTree*)f->Get("esdTree");
cd259404	67	AliESDEvent* esd = new AliESDEvent();
cd259404	68	esd->ReadFromTree(esdTree);
77e93dc2	69
c8fe2783	70	//Init geometry and array that will contain the clusters
	71	TRefArray* caloClusters = new TRefArray();
	72	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEAR") ;
	73	Float_t pos[3];
	74
	75	//Loop of events
cd259404	76	Int_t nEvt = esdTree->GetEntries();
cd259404	77	for(Int_t iev = 0; iev < nEvt; iev++) {
9480984d	78	cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
cd259404	79	esdTree->GetEvent(iev);
c8fe2783	80
	81	//Pass the geometry transformation matrix from ESDs to geometry
	82	for(Int_t mod=0; mod < (geom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
	83	//printf("matrix %d\n",mod);
	84	if(esd->GetEMCALMatrix(mod)) {
	85	//printf("EMCAL: mod %d, matrix %p\n",mod, esd->GetEMCALMatrix(mod));
	86	//(esd->GetEMCALMatrix(mod))->Print("");
	87	geom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
	88	}//matrix
	89	}//module
	90
	91	//In case you want to play with MC data, get stack
	92	AliStack *stack = 0;
	93	if(kPrintKine){
	94	AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(), "read");
	95	rl->LoadKinematics();
	96	rl->GetEvent(iev);
	97	stack = rl->Stack();
	98	}
	99
	100	//get reconstructed vertex position
	101	Double_t vertex_position[3];
	102	esd->GetVertex()->GetXYZ(vertex_position);
	103
	104	//------------------------------------------------------
	105	//Get Cells Array, all cells in event, print cells info
	106	//------------------------------------------------------
	107	AliVCaloCells &cells= *(esd->GetEMCALCells());
77e93dc2	108
c8fe2783	109	if(kPrintCaloCells){
	110	Int_t nTotalCells = cells.GetNumberOfCells() ;
	111	//Int_t type = cells.GetType();
	112	for (Int_t icell= 0; icell < nTotalCells; icell++) {
77e93dc2	113	cout<<"Cell : "<<icell<<"/"<<nTotalCells<<" - ID: "<<cells.GetCellNumber(icell)<<"; Amplitude: "<<cells.GetAmplitude(icell)<<"; Time: "<<cells.GetTime(icell)*1e9;
	114	cout << "; MC label "<<cells.GetMCLabel(icell)<<"; Embeded E fraction "<<cells.GetEFraction(icell);
	115	cout<<endl; }// cell loop
c8fe2783	116	}
c8fe2783	117
754a26d5	118	//------------------------------------------------------
	119	// Calo Trigger
	120	//------------------------------------------------------
77e93dc2	121
754a26d5	122	if(kPrintCaloTrigger)
	123	{
	124	AliESDCaloTrigger& trg = *(esd->GetCaloTrigger("EMCAL"));
	125
	126	trg.Reset();
	127	while (trg.Next())
	128	{
	129	int posX, posY;
	130	trg.GetPosition(posX, posY);
	131
	132	if (posX > -1 && posY > -1)
	133	{
	134	Int_t ts = 0;
	135	trg.GetL1TimeSum(ts);
	136
	137	cout << "Position: " << posX << " " << posY << " L1 amplitude: " << ts << endl;
	138	}
	139	}
	140	}
77e93dc2	141
c8fe2783	142	//------------------------------------------------------
	143	// Calo Clusters
	144	//------------------------------------------------------
	145
	146	//Get CaloClusters Array
c8fe2783	147	esd->GetEMCALClusters(caloClusters);
c8fe2783	148
cd259404	149	//loop over clusters
	150	Int_t nclus = caloClusters->GetEntries();
	151	for (Int_t icl = 0; icl < nclus; icl++) {
c8fe2783	152
c8fe2783	153	AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
cd259404	154	Float_t energy = clus->E();
	155	clus->GetPosition(pos);
	156	TVector3 vpos(pos[0],pos[1],pos[2]);
c8fe2783	157
	158	//We can get a momentum TLorentzVector per cluster, corrected by the vertex position
	159	//TLorentzVector p;
	160	//clus->GetMomentum(p,vertex_position);
	161
cd259404	162	Double_t cphi = vpos.Phi();
d7f5c01a	163	if(cphi < 0) cphi +=TMath::TwoPi();
cd259404	164	Double_t ceta = vpos.Eta();
c8fe2783	165
	166	Int_t nMatched = clus->GetNTracksMatched();
	167	Int_t trackIndex = clus->GetTrackMatchedIndex();
	168	Int_t nLabels = clus->GetNLabels();
cd259404	169	Int_t labelIndex = clus->GetLabel();
c8fe2783	170	Int_t nCells = clus->GetNCells();
77e93dc2	171
c8fe2783	172	//Fill some histograms
	173	hEta->Fill(ceta);
	174	hPhi->Fill(cphi*TMath::RadToDeg());
	175	hE ->Fill(energy);
	176	hTime->Fill(clus->GetTOF()*1e9);
	177
	178	//Print basic cluster information
	179	cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << "; Phi: "
77e93dc2	180	<< cphi*TMath::RadToDeg() << "; Eta: " << ceta << "; NCells: " << nCells << " ;#Matches: " << nMatched
	181	<< "; Index: " << trackIndex << "; #Labels: " << nLabels << " Index: "
	182	<< labelIndex << "; Time "<<clus->GetTOF()*1e9<<" ns "<<endl;
c8fe2783	183
	184	//Print primary info
	185	if(stack && kPrintKine) {
77e93dc2	186	if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
	187	TParticle * particle = stack->Particle(labelIndex);
	188	//Fill histograms with primary info
	189	hMCEta->Fill(particle->Eta());
	190	hMCPhi->Fill(particle->Phi()*TMath::RadToDeg());
	191	hMCE ->Fill(particle->Energy());
	192	hDEta ->Fill(ceta-particle->Eta());
	193	hDPhi ->Fill(cphiTMath::RadToDeg()-particle->Phi()TMath::RadToDeg());
	194	hDE ->Fill(energy-particle->Energy());
	195	//Print primary values
	196	cout<<" More contributing primary: "<<particle->GetName()<<"; with kinematics: "<<endl;
	197	cout<<" \t Energy: "<<particle->Energy()<<"; Phi: "<<particle->Phi()*TMath::RadToDeg()<<"; Eta: "<<particle->Eta()<<endl;
	198	for(Int_t i = 1; i < nLabels; i++){
	199	//particle = stack->Particle((((AliESDCaloCluster*)clus)->GetLabelsArray())->At(i));
	200	particle = stack->Particle((clus->GetLabels())[i]);
	201	//or Int_t *labels = clus->GetLabels();
	202	//particle = stack->Particle(labels[i]);
	203	cout<<" Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
	204	}
	205	}
	206	else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
	207	else cout<<"Negative label!!! : "<<labelIndex<<endl;
c8fe2783	208	} // play with stack
	209
	210	// Matching results
	211	if(kPrintTrackMatches && trackIndex >= 0) {
77e93dc2	212	AliESDtrack* track = esd->GetTrack(trackIndex);
	213	Double_t tphi = track->GetOuterParam()->Phi();
	214	Double_t teta = track->GetOuterParam()->Eta();
	215	Double_t tmom = track->GetOuterParam()->P();
	216	cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
	217
	218	Double_t deta = teta - ceta;
	219	Double_t dphi = tphi - cphi;
	220	if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
	221	if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
	222	Double_t dR = sqrt(dphidphi + detadeta);
	223
	224	Double_t pOverE = tmom/energy;
	225
	226	if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
	227	cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
	228	}
c8fe2783	229	}// matching
	230
	231	//Get PID weights and print them
	232	if(kPrintClusterPID){
77e93dc2	233	const Double_t *pid = clus->GetPID();
	234	printf("PID weights: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
	235	pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
	236	pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
	237	pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
	238	pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
c8fe2783	239	}//PID
	240
	241	//Get CaloCells of cluster and print their info, position.
	242	if(kPrintClusterCells){
77e93dc2	243	UShort_t * index = clus->GetCellsAbsId() ;
	244	Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
	245	Int_t sm = -1;
	246	for(Int_t i = 0; i < nCells ; i++){
	247	Int_t absId = index[i]; // or clus->GetCellNumber(i) ;
	248	Double_t ampFract = fraction[i];
	249	Float_t amp = cells.GetCellAmplitude(absId) ;
	250	Double_t time = cells.GetCellTime(absId);
	251	cout<<"\t Cluster Cell: AbsID : "<< absId << " == "<<clus->GetCellAbsId(i) <<"; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time*1e9<<endl;
	252	//Geometry methods
	253	Int_t iSupMod = 0 ;
	254	Int_t iTower = 0 ;
	255	Int_t iIphi = 0 ;
	256	Int_t iIeta = 0 ;
	257	Int_t iphi = 0 ;
	258	Int_t ieta = 0 ;
	259	if(geom){
	260	geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
	261	//Gives SuperModule and Tower numbers
	262	geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
	263	iIphi, iIeta,iphi,ieta);
	264	//Gives label of cell in eta-phi position per each supermodule
	265	Float_t cellPhi = 0;
	266	Float_t cellEta = 0;
	267	geom->EtaPhiFromIndex(absId,cellEta,cellPhi);
	268	cout<< " SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
	269	<<"; Phi "<<iIphi<<"; Index: Cell Eta "<<ieta<<"; Cell Phi "<<iphi
	270	<<"; Global: Cell Eta "<<cellEta<<"; Cell Phi "<<cellPhi*TMath::RadToDeg()<<endl;
	271	if(i==0) sm = iSupMod;
	272	else{
	273	if(sm!=iSupMod) printf("******CLUSTER SHARED BY 2 SuperModules!!!!\n");
	274	}
	275	}// geometry on
	276	}// cluster cell loop
c8fe2783	277	}// print cell clusters
	278	} //cluster loop
	279	} // event loop
	280
	281
	282	//Write histograms in a file
	283	TFile * fhisto = (TFile*) new TFile("histos.root","recreate");
	284	hEta->Write();
	285	hPhi->Write();
	286	hE ->Write();
	287	hTime->Write();
	288	if(kPrintKine){
	289	hMCEta->Write();
	290	hMCPhi->Write();
	291	hMCE ->Write();
	292	hDEta->Write();
	293	hDPhi->Write();
	294	hDE ->Write();
cd259404	295	}
c8fe2783	296	fhisto->Close();
cd259404	297	}