EMCAL/macros/TestESD.C

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