Adding AliAnalysisTaskEMCALIsolation to PWGGA/EMCALTasks
authormmarquar <marco.marquard@cern.ch>
Fri, 28 Mar 2014 16:38:23 +0000 (17:38 +0100)
committerloizides <cloizides@lbl.gov>
Fri, 28 Mar 2014 16:52:10 +0000 (17:52 +0100)
Signed-off-by: loizides <cloizides@lbl.gov>

PWGGA/CMakelibPWGGAEMCALTasks.pkg
PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx [new file with mode: 0644]
PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.h [new file with mode: 0644]
PWGGA/PWGGAEMCALTasksLinkDef.h

index 62c3410..6297131 100644 (file)
@@ -41,6 +41,7 @@ set ( SRCS
  EMCALTasks/AliAnalysisTaskEMCALCaloTrackCorr.cxx
  EMCALTasks/AliAnalysisTaskEMCALPi0V2ShSh.cxx
  EMCALTasks/AliEMCalpi0ClusterEvaluationTask.cxx
+ EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx
new file mode 100644 (file)
index 0000000..c268183
--- /dev/null
@@ -0,0 +1,489 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TObjArray.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+
+#include "AliVEvent.h"
+#include "AliVCaloCells.h"
+#include "AliVCluster.h"
+
+#include "AliEMCALGeometry.h"
+#include "AliGeomManager.h"
+
+#include "AliAnalysisTaskEMCALIsolation.h"
+
+#include "AliStack.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+
+#include "TDatabasePDG.h"
+#include "TParticlePDG.h"
+
+// Task for isolating gammas with EMCAL
+// Author: Marco Marquard
+
+ClassImp(AliAnalysisTaskEMCALIsolation)
+
+       //________________________________________________________________________
+
+AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation()
+       : AliAnalysisTaskSE(),
+       bVerbose(0),
+       bMC(0),
+       fESD(0),
+       fAOD(0),
+       fMC(0),
+       fStack(0),
+       fOutputList(0),
+       fHistGlobalHmap(0),
+       fHistGlobalHmap0(0),
+       fTreeEvent(0),
+       fTreeCluster(0),
+       emclus(0),
+       prodrad(0),
+       contPID(-999),
+       mothPID(-999),
+       trackmatch(0),
+       fGeom(0),
+       fGeoName("EMCAL_COMPLETE"),
+       fESDCells(0),
+       fAODCells(0),
+       fEsdClusters(0),
+       fAodClusters(0)
+{
+       //default constructor
+       DefineInput(0, TChain::Class());
+       DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation(const char *name) 
+       : AliAnalysisTaskSE(name),
+       bVerbose(0),
+       bMC(0),
+       fESD(0),
+       fAOD(0),
+       fMC(0),
+       fStack(0),
+       fOutputList(0),
+       fHistGlobalHmap(0),
+       fHistGlobalHmap0(0),
+       fTreeEvent(0),
+       fTreeCluster(0),
+       emclus(0),
+       prodrad(0),
+       contPID(-999),
+       mothPID(-999),
+       trackmatch(0),
+       fGeom(0),
+       fGeoName("EMCAL_COMPLETE"),
+       fESDCells(0),
+       fAODCells(0),
+       fEsdClusters(0),
+       fAodClusters(0)
+{
+       // Constructor
+
+       // Define input and output slots here
+       // Input slot #0 works with a TChain
+       DefineInput(0, TChain::Class());
+       // Output slot #0 id reserved by the base class for AOD
+       // Output slot #1 writes into a TH1 container
+       DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEMCALIsolation::~AliAnalysisTaskEMCALIsolation() 
+{
+       //Destructor
+
+       delete fESD;
+       delete fAOD;
+       delete fOutputList;
+       delete fHistGlobalHmap;
+       delete fHistGlobalHmap0;
+       delete fTreeCluster;
+       delete fTreeEvent;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALIsolation::UserCreateOutputObjects()
+{
+       // Create histograms
+       // Called once
+       Double_t etagaps[97]= {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
+
+       Double_t phigaps[125]= {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732,1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
+
+       fOutputList = new TList();
+
+       fHistGlobalHmap = new TH2F("fHistGlobalHmap","Hit map for EMCAL",96,etagaps,124,phigaps);
+       fHistGlobalHmap->GetXaxis()->SetTitle("#eta");
+       fHistGlobalHmap->GetYaxis()->SetTitle("#phi");
+       fOutputList->Add(fHistGlobalHmap);
+
+       fHistGlobalHmap0 = new TH2F("fHistGlobalHmap0","Hit map for EMCAL",96,0,96,120,0,120);
+       fHistGlobalHmap0->GetXaxis()->SetTitle("#eta");
+       fHistGlobalHmap0->GetYaxis()->SetTitle("#phi");
+       fOutputList->Add(fHistGlobalHmap0);
+       /*      
+               fHist = new TH1F("","",,,,);
+               fHist->GetXaxis()->SetTitle("");
+               fHist->GetYaxis()->SetTitle("");
+               fHist->SetMarkerStyle(kFullCircle);
+               fOutputList->Add(fHist);
+               */
+
+       fTreeEvent = new TTree("eventsT","eventsT");
+       fTreeEvent->Branch("emclus",    &emclus);
+
+       fTreeCluster = new TTree("clusterT","clusterT");
+       fTreeCluster->Branch("emclus",  &emclus);
+       fTreeCluster->Branch("prodrad", &prodrad);
+       fTreeCluster->Branch("contPID", &contPID);
+       fTreeCluster->Branch("mothPID", &mothPID);
+
+       fOutputList->Add(fHistGlobalHmap);
+       fOutputList->Add(fHistGlobalHmap0);
+       fOutputList->Add(fTreeEvent);
+       fOutputList->Add(fTreeCluster);
+
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALIsolation::UserExec(Option_t *) 
+{
+
+       AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+
+
+       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
+       if (!esdH)
+       {
+               Printf("ERROR: Could not get ESDInputHandler");
+       }
+
+       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
+       if (!aodH && !esdH)
+       {
+               Printf("ERROR: Could not get AODInputHandler");
+       }
+
+       if(esdH){
+               fESD = esdH->GetEvent();
+               am->LoadBranch("AliESDRun.");
+               am->LoadBranch("AliESDHeader.");
+       }
+       else if(aodH){
+               fAOD = aodH->GetEvent();
+       }
+       else{
+               AliFatal("Neither ESD nor AOD event found");
+               return;
+       }
+
+       AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+       if (!mcH) {
+               Printf("ERROR: Could not retrieve MC event handler");
+               return;
+       }
+       fMC = mcH->MCEvent();
+       if (!fMC) {
+               Printf("ERROR: Could not retrieve MC event");
+               return;
+       }
+       // stack to get truth infos of particles:
+       fStack = fMC->Stack();
+       if (!fStack) {
+               Printf("ERROR: Could not retrieve stack from MC event");
+               return;
+       }
+       //Printf(" MC particles: %d", fMC->GetNumberOfTracks());
+
+
+       fGeom = AliEMCALGeometry::GetInstance(fGeoName);
+       //AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
+       //Double_t phimin = emc->GetArm1PhiMin();
+       //Double_t phimax = emc->GetArm1PhiMax();
+       //printf("EMCAL coverage in phi from: %f to: %f \n",phimin,phimax);
+       if (!fGeom){
+               AliFatal("EMCAL(fGeom) geometry is not initialized");
+               return;
+       }
+
+       AliVCaloCells *cells = 0;
+
+       if (fESD){
+               am->LoadBranch("EMCALCells."); 
+               fESDCells = fESD->GetEMCALCells();
+               cells = fESDCells;
+       }
+       else{
+               am->LoadBranch("emcalCells");
+               fAODCells = fAOD->GetEMCALCells();
+               cells = fAODCells;
+       }
+
+
+
+
+       if (!cells){
+               AliFatal("No cells loaded");
+               return;
+       }
+
+       if(fESD){ 
+               am->LoadBranch("CaloClusters");
+               TList *l = fESD->GetList();
+               fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
+       }
+       else if(fAOD){
+               am->LoadBranch("caloClusters");
+               TList *l = fAOD->GetList();
+               fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
+       }
+       else{
+               AliFatal("Impossible to not have either pointer to ESD or AOD event");
+       }
+
+       TObjArray *clusters = fEsdClusters;
+       if (!clusters)
+               clusters = fAodClusters;
+       if (!clusters){
+               Printf("ERROR: no clusters node in event!");
+               //return;
+       }
+
+       //========================================= start of analysis ============================================
+
+
+       //============> DATA LOOP =================>
+
+       //------------> cell loop ----------------->
+       Int_t ncells = cells->GetNumberOfCells();
+       //for (Int_t i = 0; i<15000; ++i){
+       //      Int_t absID    = i;
+       for (Int_t i = 0; i<ncells; ++i){
+               Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
+
+               //Check if this absId exists
+               if(!(fGeom->CheckAbsCellId(absID))){ 
+                       continue;
+               }
+               Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
+               // Get from the absId the super module number, the module number and the eta-phi index (0 or 1) in the module
+               fGeom->GetCellIndex(absID, nSupMod, nModule, nIphi, nIeta);
+
+               // Get from the  the super module number, the module number and the eta-phi index (0 or 1) in the module the tower row (iphi) and column (ieta)
+               fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
+
+               if(bVerbose)
+                       printf("ID: %d ; supermodule: %d ; module %d ; phi %d ; eta %d \n",absID, nSupMod, nModule, nIphi, nIeta);
+
+
+               Int_t etainv = 47-ieta;
+               //fHSMCellE[nSupMod]->Fill(etainv,phi,i);
+
+
+               Int_t binx = etainv + 48*((nSupMod+1) % 2);
+               Int_t biny = iphi + 24*nSupMod/2;
+               if(nSupMod % 2 != 0)
+                       biny -= 12;
+
+               Double_t cellE = cells->GetAmplitude(i);
+               //fHistGlobalHmap0->Fill(binx,biny,cellE);
+
+               Double_t nEtaGlobal, nPhiGlobal; 
+               Int_t nSupMod2;
+
+               nSupMod2 = fGeom->GetSuperModuleNumber(absID);
+               fGeom->EtaPhiFromIndex(absID,nEtaGlobal,nPhiGlobal);
+
+               if(bVerbose){
+                       printf("SM number: %d ; global eta: \t%6.4f ; global phi: \t%6.4f \n",nSupMod2,nEtaGlobal,nPhiGlobal);
+               }
+
+
+               //Double_t nPhi;
+               //if(nPhiGlobal < 0){
+               //      nPhi = TMath::Pi()+TMath::Abs(TMath::Pi()+nPhiGlobal);
+               //}
+               //else{
+               //nPhi = nPhiGlobal;
+               //}
+               fHistGlobalHmap0->Fill(binx,biny,cellE);
+               fHistGlobalHmap->Fill(nEtaGlobal,nPhiGlobal,cellE);
+
+       }
+
+       if(bVerbose){
+               printf("\nThere are %d cells in this event\n", ncells);
+               printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());
+       }
+
+
+       //<------------ cell loop <-----------------
+
+       //------------> track loop ----------------->
+       for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
+               AliESDtrack* track = fESD->GetTrack(iTracks);
+               if (!track) {
+                       printf("ERROR: Could not receive track %d\n", iTracks);
+                       continue;
+               }
+
+
+       }
+       //<------------ track loop <-----------------
+
+       //------------> cluster loop ----------------->
+       emclus = 0;
+       Int_t nclus = clusters->GetEntries();
+       //printf("%d clusters\n",nclus);
+       //cluster loop 1 for counting EMCal cluster
+       for(int i = 0; i < nclus; i++){ 
+               AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
+               if (!clus){
+                       //printf("Cluster empty!\n");
+                       continue;
+               }
+               if (!clus->IsEMCAL()){
+                       //printf("Cluster %d is not EMCAL!\n",i);
+                       //      cout << "at " << i << " cluster not EMCAL" << endl;
+                       continue;
+               }
+
+               emclus++;
+       }
+
+       //<------------ cluster loop <-----------------
+
+
+       //<============ DATA LOOP <=================
+
+
+       //============> MC LOOP =================>
+       if(bMC){
+
+               //------------> cell loop ----------------->
+               //<------------ cell loop <-----------------
+
+               //------------> track loop ----------------->
+               //<------------ track loop <-----------------
+
+               //------------> cluster loop  ----------------->
+               for(int i = 0; i < nclus; i++){
+                       AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
+                       if (!clus){
+                               //printf("Cluster empty!\n");
+                               continue;
+                       }
+                       if (!clus->IsEMCAL()){
+                               //printf("Cluster %d is not EMCAL!\n",i);
+                               //      cout << "at " << i << " cluster not EMCAL" << endl;
+                               continue;
+                       }
+                       //MCTruth
+                       int clslabel = -1;
+
+                       // all contributors
+                       Int_t* mcarr =  clus->GetLabels();
+                       // number of contributors
+                       Int_t nl = clus->GetNLabels();
+
+
+
+
+                       for(int ip=0;ip<nl;ip++){
+                               Int_t entry = mcarr[ip];
+                               AliMCParticle *mcPart = static_cast<AliMCParticle*>(fMC->GetTrack(entry));
+                               if(!mcPart){
+                                       continue;
+                               }
+                               contPID = mcPart->PdgCode();
+                               if(bVerbose){
+                                       printf("%d \t %d \n",clslabel,mcPart->PdgCode());
+                               }
+                               prodrad = mcPart->Particle()->R();
+                               Int_t indexMoth = mcPart->GetMother();
+                               if(indexMoth >= 0)
+                               {
+                                       AliMCParticle* moth = static_cast<AliMCParticle*>(fMC->GetTrack(indexMoth));
+                                       //      //AliMCParticle* moth = (AliMCParticle*)stack->At(indexMoth);
+                                       mothPID = moth->PdgCode();
+                               }
+
+
+
+
+
+                               //printf("%d \t %d \n",clslabel,nl);
+
+                       }
+                       fTreeCluster->Fill();
+
+               }
+
+               //<------------ cluster loop <-----------------
+
+       }
+       //<============ MC LOOP <=================
+
+
+
+
+
+
+
+       //cluster loop 2 for analysis
+       //\cluster loop
+
+
+       fTreeEvent->Fill();
+
+
+       PostData(1, fOutputList);
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALIsolation::Terminate(Option_t *) 
+{
+       // Draw result to the screen
+       // Called once at the end of the query
+
+       fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+       if (!fOutputList) {
+               printf("ERROR: Output list not available\n");
+               return;
+       }
+
+       fHistGlobalHmap = dynamic_cast<TH2F*> (fOutputList->At(2));
+       if (!fHistGlobalHmap) {
+               printf("ERROR: fHistGlobalHmap not available\n");
+               return;
+       }
+
+       fHistGlobalHmap->DrawCopy("colz");
+
+}
+
+
+const char * AliAnalysisTaskEMCALIsolation::GetParticleName(Int_t pdg)
+{
+       TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
+       if(p1) return p1->GetName();
+       return Form("%d", pdg);
+}
diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.h b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.h
new file mode 100644 (file)
index 0000000..6b58922
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef AliAnalysisTaskEMCALIsolation_cxx
+#define AliAnalysisTaskEMCALIsolation_cxx
+
+// Task for isolating gammas with EMCAL
+// Author: Marco Marquard
+
+class TH1F;
+class AliESDEvent;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskEMCALIsolation : public AliAnalysisTaskSE {
+       public:
+               AliAnalysisTaskEMCALIsolation();
+               AliAnalysisTaskEMCALIsolation(const char *name);
+               virtual ~AliAnalysisTaskEMCALIsolation();
+
+               virtual void   UserCreateOutputObjects();
+               virtual void   UserExec(Option_t *option);
+               virtual void   Terminate(Option_t *);
+
+               void    SetVerbose(Bool_t b)                    {bVerbose       = b;}
+               void    SetMC(Bool_t mc)                        {bMC    = mc;}
+
+       protected:
+
+               const char * GetParticleName(Int_t);    
+
+       private:
+               Bool_t          bVerbose;                       // Verbose option
+               Bool_t          bMC;                            // MC option
+               AliESDEvent     *fESD;                          // ESD object
+               AliAODEvent     *fAOD;                          // AOD object
+               AliMCEvent      *fMC;                           // MC Event
+               AliStack        *fStack;                        // Ali stack
+               TList           *fOutputList;                   // Output list
+               TH2F            *fHistGlobalHmap;               // Cell hit map for the complete EMCAL
+               TH2F            *fHistGlobalHmap0;              // Cell hit map for the complete EMCAL
+               TTree           *fTreeEvent;                    //Tree with event informations
+               TTree           *fTreeCluster;                  //Tree with cluster informations
+               Int_t           emclus;                         //number of cluster per event
+               Double_t        prodrad;                        //production radius of V0 vertex
+               Int_t           contPID;                        //PID of contributor
+               Int_t           mothPID;                        //PID of contributor mother
+               Bool_t          trackmatch;                             //Track matching
+
+               AliEMCALGeometry        *fGeom;                 // geometry utils
+               TString                 fGeoName;               // geometry name (def = EMCAL_COMPLET (alternative: EMCAL_FIRSTYEARV1))
+               AliESDCaloCells *fESDCells;                     //!pointer to esd cells
+               AliAODCaloCells *fAODCells;                     //!pointer to aod cells
+
+               TObjArray       *fEsdClusters;                  //!pointer to esd clusters
+               TObjArray       *fAodClusters;                  //!pointer to aod clusters
+
+
+               AliAnalysisTaskEMCALIsolation(const AliAnalysisTaskEMCALIsolation&); // not implemented
+               AliAnalysisTaskEMCALIsolation& operator=(const AliAnalysisTaskEMCALIsolation&); // not implemented
+
+               ClassDef(AliAnalysisTaskEMCALIsolation, 1); // example of analysis
+};
+
+#endif
index 7db82e4..82f29e4 100644 (file)
@@ -31,5 +31,6 @@
 #pragma link C++ class AliAnalysisTaskEMCALCaloTrackCorr+;
 #pragma link C++ class AliAnalysisTaskEMCALPi0V2ShSh+;
 #pragma link C++ class AliEMCalpi0ClusterEvaluationTask+;
+#pragma link C++ class AliAnalysisTaskEMCALIsolation+;
 
 #endif