PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TCanvas.h"
   5 #include "TFile.h"
   6 #include "TObjArray.h"
   7
   8 #include "AliAnalysisTask.h"
   9 #include "AliAnalysisManager.h"
  10
  11 #include "AliESDEvent.h"
  12 #include "AliESDInputHandler.h"
  13 #include "AliAODEvent.h"
  14 #include "AliAODInputHandler.h"
  15
  16 #include "AliVEvent.h"
  17 #include "AliVCaloCells.h"
  18 #include "AliVCluster.h"
  19
  20 #include "AliEMCALGeometry.h"
  21 #include "AliGeomManager.h"
  22
  23 #include "AliAnalysisTaskEMCALIsolation.h"
  24
  25 #include "AliStack.h"
  26 #include "AliMCEventHandler.h"
  27 #include "AliMCEvent.h"
  28 #include "AliMCParticle.h"
  29
  30 #include "TDatabasePDG.h"
  31 #include "TParticlePDG.h"
  32
  33 // Task for isolating gammas with EMCAL
  34 // Author: Marco Marquard
  35
  36 ClassImp(AliAnalysisTaskEMCALIsolation)
  37
  38         //________________________________________________________________________
  39
  40 AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation()
  41         : AliAnalysisTaskSE(),
  42         bVerbose(0),
  43         bMC(0),
  44         fESD(0),
  45         fAOD(0),
  46         fMC(0),
  47         fStack(0),
  48         fOutputList(0),
  49         fHistGlobalHmap(0),
  50         fHistGlobalHmap0(0),
  51         fTreeEvent(0),
  52         fTreeCluster(0),
  53         emclus(0),
  54         prodrad(0),
  55         contPID(-999),
  56         mothPID(-999),
  57         trackmatch(0),
  58         fGeom(0),
  59         fGeoName("EMCAL_COMPLETE"),
  60         fESDCells(0),
  61         fAODCells(0),
  62         fEsdClusters(0),
  63         fAodClusters(0)
  64 {
  65         //default constructor
  66         DefineInput(0, TChain::Class());
  67         DefineOutput(1, TList::Class());
  68 }
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation(const char *name)
  72         : AliAnalysisTaskSE(name),
  73         bVerbose(0),
  74         bMC(0),
  75         fESD(0),
  76         fAOD(0),
  77         fMC(0),
  78         fStack(0),
  79         fOutputList(0),
  80         fHistGlobalHmap(0),
  81         fHistGlobalHmap0(0),
  82         fTreeEvent(0),
  83         fTreeCluster(0),
  84         emclus(0),
  85         prodrad(0),
  86         contPID(-999),
  87         mothPID(-999),
  88         trackmatch(0),
  89         fGeom(0),
  90         fGeoName("EMCAL_COMPLETE"),
  91         fESDCells(0),
  92         fAODCells(0),
  93         fEsdClusters(0),
  94         fAodClusters(0)
  95 {
  96         // Constructor
  97
  98         // Define input and output slots here
  99         // Input slot #0 works with a TChain
 100         DefineInput(0, TChain::Class());
 101         // Output slot #0 id reserved by the base class for AOD
 102         // Output slot #1 writes into a TH1 container
 103         DefineOutput(1, TList::Class());
 104 }
 105
 106 //________________________________________________________________________
 107 AliAnalysisTaskEMCALIsolation::~AliAnalysisTaskEMCALIsolation()
 108 {
 109         //Destructor
 110
 111         delete fESD;
 112         delete fAOD;
 113         delete fOutputList;
 114         delete fHistGlobalHmap;
 115         delete fHistGlobalHmap0;
 116         delete fTreeCluster;
 117         delete fTreeEvent;
 118 }
 119
 120 //________________________________________________________________________
 121 void AliAnalysisTaskEMCALIsolation::UserCreateOutputObjects()
 122 {
 123         // Create histograms
 124         // Called once
 125         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};
 126
 127         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};
 128
 129         fOutputList = new TList();
 130
 131         fHistGlobalHmap = new TH2F("fHistGlobalHmap","Hit map for EMCAL",96,etagaps,124,phigaps);
 132         fHistGlobalHmap->GetXaxis()->SetTitle("#eta");
 133         fHistGlobalHmap->GetYaxis()->SetTitle("#phi");
 134         fOutputList->Add(fHistGlobalHmap);
 135
 136         fHistGlobalHmap0 = new TH2F("fHistGlobalHmap0","Hit map for EMCAL",96,0,96,120,0,120);
 137         fHistGlobalHmap0->GetXaxis()->SetTitle("#eta");
 138         fHistGlobalHmap0->GetYaxis()->SetTitle("#phi");
 139         fOutputList->Add(fHistGlobalHmap0);
 140         /*
 141                 fHist = new TH1F("","",,,,);
 142                 fHist->GetXaxis()->SetTitle("");
 143                 fHist->GetYaxis()->SetTitle("");
 144                 fHist->SetMarkerStyle(kFullCircle);
 145                 fOutputList->Add(fHist);
 146                 */
 147
 148         fTreeEvent = new TTree("eventsT","eventsT");
 149         fTreeEvent->Branch("emclus",    &emclus);
 150
 151         fTreeCluster = new TTree("clusterT","clusterT");
 152         fTreeCluster->Branch("emclus",  &emclus);
 153         fTreeCluster->Branch("prodrad", &prodrad);
 154         fTreeCluster->Branch("contPID", &contPID);
 155         fTreeCluster->Branch("mothPID", &mothPID);
 156
 157         fOutputList->Add(fHistGlobalHmap);
 158         fOutputList->Add(fHistGlobalHmap0);
 159         fOutputList->Add(fTreeEvent);
 160         fOutputList->Add(fTreeCluster);
 161
 162
 163
 164 }
 165
 166 //________________________________________________________________________
 167 void AliAnalysisTaskEMCALIsolation::UserExec(Option_t *)
 168 {
 169
 170         AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 171
 172
 173         AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
 174         if (!esdH)
 175         {
 176                 Printf("ERROR: Could not get ESDInputHandler");
 177         }
 178
 179         AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
 180         if (!aodH && !esdH)
 181         {
 182                 Printf("ERROR: Could not get AODInputHandler");
 183         }
 184
 185         if(esdH){
 186                 fESD = esdH->GetEvent();
 187                 am->LoadBranch("AliESDRun.");
 188                 am->LoadBranch("AliESDHeader.");
 189         }
 190         else if(aodH){
 191                 fAOD = aodH->GetEvent();
 192         }
 193         else{
 194                 AliFatal("Neither ESD nor AOD event found");
 195                 return;
 196         }
 197
 198         AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 199         if (!mcH) {
 200                 Printf("ERROR: Could not retrieve MC event handler");
 201                 return;
 202         }
 203         fMC = mcH->MCEvent();
 204         if (!fMC) {
 205                 Printf("ERROR: Could not retrieve MC event");
 206                 return;
 207         }
 208         // stack to get truth infos of particles:
 209         fStack = fMC->Stack();
 210         if (!fStack) {
 211                 Printf("ERROR: Could not retrieve stack from MC event");
 212                 return;
 213         }
 214         //Printf(" MC particles: %d", fMC->GetNumberOfTracks());
 215
 216
 217         fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 218         //AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 219         //Double_t phimin = emc->GetArm1PhiMin();
 220         //Double_t phimax = emc->GetArm1PhiMax();
 221         //printf("EMCAL coverage in phi from: %f to: %f \n",phimin,phimax);
 222         if (!fGeom){
 223                 AliFatal("EMCAL(fGeom) geometry is not initialized");
 224                 return;
 225         }
 226
 227         AliVCaloCells *cells = 0;
 228
 229         if (fESD){
 230                 am->LoadBranch("EMCALCells.");
 231                 fESDCells = fESD->GetEMCALCells();
 232                 cells = fESDCells;
 233         }
 234         else{
 235                 am->LoadBranch("emcalCells");
 236                 fAODCells = fAOD->GetEMCALCells();
 237                 cells = fAODCells;
 238         }
 239
 240
 241
 242
 243         if (!cells){
 244                 AliFatal("No cells loaded");
 245                 return;
 246         }
 247
 248         if(fESD){
 249                 am->LoadBranch("CaloClusters");
 250                 TList *l = fESD->GetList();
 251                 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 252         }
 253         else if(fAOD){
 254                 am->LoadBranch("caloClusters");
 255                 TList *l = fAOD->GetList();
 256                 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 257         }
 258         else{
 259                 AliFatal("Impossible to not have either pointer to ESD or AOD event");
 260         }
 261
 262         TObjArray *clusters = fEsdClusters;
 263         if (!clusters)
 264                 clusters = fAodClusters;
 265         if (!clusters){
 266                 Printf("ERROR: no clusters node in event!");
 267                 //return;
 268         }
 269
 270         //========================================= start of analysis ============================================
 271
 272
 273         //============> DATA LOOP =================>
 274
 275         //------------> cell loop ----------------->
 276         Int_t ncells = cells->GetNumberOfCells();
 277         //for (Int_t i = 0; i<15000; ++i){
 278         //      Int_t absID    = i;
 279         for (Int_t i = 0; i<ncells; ++i){
 280                 Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
 281
 282                 //Check if this absId exists
 283                 if(!(fGeom->CheckAbsCellId(absID))){
 284                         continue;
 285                 }
 286                 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
 287                 // Get from the absId the super module number, the module number and the eta-phi index (0 or 1) in the module
 288                 fGeom->GetCellIndex(absID, nSupMod, nModule, nIphi, nIeta);
 289
 290                 // 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)
 291                 fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
 292
 293                 if(bVerbose)
 294                         printf("ID: %d ; supermodule: %d ; module %d ; phi %d ; eta %d \n",absID, nSupMod, nModule, nIphi, nIeta);
 295
 296
 297                 Int_t etainv = 47-ieta;
 298                 //fHSMCellE[nSupMod]->Fill(etainv,phi,i);
 299
 300
 301                 Int_t binx = etainv + 48*((nSupMod+1) % 2);
 302                 Int_t biny = iphi + 24*nSupMod/2;
 303                 if(nSupMod % 2 != 0)
 304                         biny -= 12;
 305
 306                 Double_t cellE = cells->GetAmplitude(i);
 307                 //fHistGlobalHmap0->Fill(binx,biny,cellE);
 308
 309                 Double_t nEtaGlobal, nPhiGlobal;
 310                 Int_t nSupMod2;
 311
 312                 nSupMod2 = fGeom->GetSuperModuleNumber(absID);
 313                 fGeom->EtaPhiFromIndex(absID,nEtaGlobal,nPhiGlobal);
 314
 315                 if(bVerbose){
 316                         printf("SM number: %d ; global eta: \t%6.4f ; global phi: \t%6.4f \n",nSupMod2,nEtaGlobal,nPhiGlobal);
 317                 }
 318
 319
 320                 //Double_t nPhi;
 321                 //if(nPhiGlobal < 0){
 322                 //      nPhi = TMath::Pi()+TMath::Abs(TMath::Pi()+nPhiGlobal);
 323                 //}
 324                 //else{
 325                 //nPhi = nPhiGlobal;
 326                 //}
 327                 fHistGlobalHmap0->Fill(binx,biny,cellE);
 328                 fHistGlobalHmap->Fill(nEtaGlobal,nPhiGlobal,cellE);
 329
 330         }
 331
 332         if(bVerbose){
 333                 printf("\nThere are %d cells in this event\n", ncells);
 334                 printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());
 335         }
 336
 337
 338         //<------------ cell loop <-----------------
 339
 340         //------------> track loop ----------------->
 341         for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 342                 AliESDtrack* track = fESD->GetTrack(iTracks);
 343                 if (!track) {
 344                         printf("ERROR: Could not receive track %d\n", iTracks);
 345                         continue;
 346                 }
 347
 348
 349         }
 350         //<------------ track loop <-----------------
 351
 352         //------------> cluster loop ----------------->
 353         emclus = 0;
 354         Int_t nclus = clusters->GetEntries();
 355         //printf("%d clusters\n",nclus);
 356         //cluster loop 1 for counting EMCal cluster
 357         for(int i = 0; i < nclus; i++){
 358                 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 359                 if (!clus){
 360                         //printf("Cluster empty!\n");
 361                         continue;
 362                 }
 363                 if (!clus->IsEMCAL()){
 364                         //printf("Cluster %d is not EMCAL!\n",i);
 365                         //      cout << "at " << i << " cluster not EMCAL" << endl;
 366                         continue;
 367                 }
 368
 369                 emclus++;
 370         }
 371
 372         //<------------ cluster loop <-----------------
 373
 374
 375         //<============ DATA LOOP <=================
 376
 377
 378         //============> MC LOOP =================>
 379         if(bMC){
 380
 381                 //------------> cell loop ----------------->
 382                 //<------------ cell loop <-----------------
 383
 384                 //------------> track loop ----------------->
 385                 //<------------ track loop <-----------------
 386
 387                 //------------> cluster loop  ----------------->
 388                 for(int i = 0; i < nclus; i++){
 389                         AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 390                         if (!clus){
 391                                 //printf("Cluster empty!\n");
 392                                 continue;
 393                         }
 394                         if (!clus->IsEMCAL()){
 395                                 //printf("Cluster %d is not EMCAL!\n",i);
 396                                 //      cout << "at " << i << " cluster not EMCAL" << endl;
 397                                 continue;
 398                         }
 399                         //MCTruth
 400                         int clslabel = -1;
 401
 402                         // all contributors
 403                         Int_t* mcarr =  clus->GetLabels();
 404                         // number of contributors
 405                         Int_t nl = clus->GetNLabels();
 406
 407
 408
 409
 410                         for(int ip=0;ip<nl;ip++){
 411                                 Int_t entry = mcarr[ip];
 412                                 AliMCParticle *mcPart = static_cast<AliMCParticle*>(fMC->GetTrack(entry));
 413                                 if(!mcPart){
 414                                         continue;
 415                                 }
 416                                 contPID = mcPart->PdgCode();
 417                                 if(bVerbose){
 418                                         printf("%d \t %d \n",clslabel,mcPart->PdgCode());
 419                                 }
 420                                 prodrad = mcPart->Particle()->R();
 421                                 Int_t indexMoth = mcPart->GetMother();
 422                                 if(indexMoth >= 0)
 423                                 {
 424                                         AliMCParticle* moth = static_cast<AliMCParticle*>(fMC->GetTrack(indexMoth));
 425                                         //      //AliMCParticle* moth = (AliMCParticle*)stack->At(indexMoth);
 426                                         mothPID = moth->PdgCode();
 427                                 }
 428
 429
 430
 431
 432
 433                                 //printf("%d \t %d \n",clslabel,nl);
 434
 435                         }
 436                         fTreeCluster->Fill();
 437
 438                 }
 439
 440                 //<------------ cluster loop <-----------------
 441
 442         }
 443         //<============ MC LOOP <=================
 444
 445
 446
 447
 448
 449
 450
 451         //cluster loop 2 for analysis
 452         //\cluster loop
 453
 454
 455         fTreeEvent->Fill();
 456
 457
 458         PostData(1, fOutputList);
 459 }
 460
 461 //________________________________________________________________________
 462 void AliAnalysisTaskEMCALIsolation::Terminate(Option_t *)
 463 {
 464         // Draw result to the screen
 465         // Called once at the end of the query
 466
 467         fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 468         if (!fOutputList) {
 469                 printf("ERROR: Output list not available\n");
 470                 return;
 471         }
 472
 473         fHistGlobalHmap = dynamic_cast<TH2F*> (fOutputList->At(2));
 474         if (!fHistGlobalHmap) {
 475                 printf("ERROR: fHistGlobalHmap not available\n");
 476                 return;
 477         }
 478
 479         fHistGlobalHmap->DrawCopy("colz");
 480
 481 }
 482
 483
 484 const char * AliAnalysisTaskEMCALIsolation::GetParticleName(Int_t pdg)
 485 {
 486         TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
 487         if(p1) return p1->GetName();
 488         return Form("%d", pdg);
 489 }