PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.cxx

   1 // $Id$
   2 //
   3 //
   4 //
   5 //
   6
   7 #include "TChain.h"
   8 #include "TTree.h"
   9 #include "TH1F.h"
  10 #include "TH2F.h"
  11 #include "TCanvas.h"
  12
  13 #include "AliAnalysisTask.h"
  14 #include "AliAnalysisManager.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliESDHeader.h"
  18 #include "AliESDUtils.h"
  19 #include "AliESDInputHandler.h"
  20 #include "AliCentrality.h"
  21 #include "AliESDpid.h"
  22 #include "AliKFParticle.h"
  23
  24 #include "AliMCEventHandler.h"
  25 #include "AliMCEvent.h"
  26 #include "AliStack.h"
  27 #include "TParticle.h"
  28 #include "AliAODMCParticle.h"
  29
  30
  31 #include "AliESDtrackCuts.h"
  32 #include "AliESDv0.h"
  33 #include "AliV0vertexer.h"
  34 #include "AliESDCaloCluster.h"
  35 #include "AliESDCaloCells.h"
  36 #include "AliAODEvent.h"
  37 #include "AliEMCALGeometry.h"
  38 #include "AliEMCALRecoUtils.h"
  39 #include "AliOADBContainer.h"
  40 #include "AliVEvent.h"
  41 #include "AliVVertex.h"
  42 #include "AliVCluster.h"
  43 #include "AliVCaloCells.h"
  44 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  45 #include "TLorentzVector.h"
  46
  47
  48 #include "AliAnalysisTaskEMCALPhoton.h"
  49 #include "TFile.h"
  50
  51
  52 using std::cout;
  53 using std::endl;
  54
  55 ClassImp(AliAnalysisTaskEMCALPhoton)
  56
  57 //________________________________________________________________________
  58 AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton() :
  59   AliAnalysisTaskSE(),
  60   fTrCuts(0),
  61   fPrTrCuts(0),
  62   fSelTracks(0),
  63   fSelPrimTracks(0),
  64   fTracks(0),
  65   fPhotConvArray(0),
  66   fMyClusts(0),
  67   fMyAltClusts(0),
  68   fMyCells(0),
  69   fMyTracks(0),
  70   fMyMcParts(0),
  71   fHeader(0x0),
  72   fOADBContainer(0),
  73   fCaloClusters(0),
  74   fCaloClustersNew(0),
  75   fAODMCParticles(0),
  76   fVCells(0),
  77   fGeom(0x0),
  78   fTimeResTOF(0),
  79   fMipResponseTPC(0),
  80   fGeoName("EMCAL_COMPLETEV1"),
  81   fPeriod("LHC11d"),
  82   fIsTrain(0),
  83   fIsMC(0),
  84   fDebug(0),
  85   fRedoV0(1),
  86   fIsGrid(0),
  87   fClusThresh(2.0),
  88   fClusterizer(0),
  89   fCaloClustersName("EmcalClusterv2"),
  90   fESD(0),
  91   fAOD(0),
  92   fVev(0),
  93   fMCEvent(0),
  94   fStack(0x0),
  95   fOutputList(0),
  96   fTree(0),
  97   fMyMcIndex(0),
  98   fNV0sBefAndAftRerun(0),
  99   fConversionVtxXY(0),
 100   fInvMassV0(0),
 101   fInvMassV0KF(0),
 102   fInvMassV0SS(0),
 103   fDedxPAll(0)
 104 {
 105   // Default constructor.
 106   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 107 }
 108
 109 //________________________________________________________________________
 110 AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton(const char *name) :
 111   AliAnalysisTaskSE(name),
 112   fTrCuts(0),
 113   fPrTrCuts(0),
 114   fSelTracks(0),
 115   fSelPrimTracks(0),
 116   fTracks(0),
 117   fPhotConvArray(0),
 118   fMyClusts(0),
 119   fMyAltClusts(0),
 120   fMyCells(0),
 121   fMyTracks(0),
 122   fMyMcParts(0),
 123   fHeader(0),
 124   fOADBContainer(0),
 125   fCaloClusters(0),
 126   fCaloClustersNew(0),
 127   fAODMCParticles(0),
 128   fVCells(0),
 129   fGeom(0x0),
 130   fTimeResTOF(0),
 131   fMipResponseTPC(0),
 132   fGeoName("EMCAL_COMPLETEV1"),
 133   fPeriod("LHC11c"),
 134   fIsTrain(0),
 135   fIsMC(0),
 136   fDebug(0),
 137   fRedoV0(1),
 138   fIsGrid(0),
 139   fClusThresh(2.),
 140   fClusterizer(0),
 141   fCaloClustersName("EmcalClusterv2"),
 142   fESD(0),
 143   fAOD(0),
 144   fVev(0),
 145   fMCEvent(0),
 146   fStack(0x0),
 147   fOutputList(0),
 148   fTree(0),
 149   fMyMcIndex(0),
 150   fNV0sBefAndAftRerun(0),
 151   fConversionVtxXY(0),
 152   fInvMassV0(0),
 153   fInvMassV0KF(0),
 154   fInvMassV0SS(0),
 155   fDedxPAll(0)
 156 {
 157   // Constructor
 158
 159   // Define input and output slots here
 160   // Input slot #0 works with a TChain
 161   DefineInput(0, TChain::Class());
 162   // Output slot #0 id reserved by the base class for AOD
 163   // Output slot #1 writes into a TH1 container
 164   DefineOutput(1, TList::Class());
 165   DefineOutput(2, TTree::Class());
 166 }
 167
 168 //________________________________________________________________________
 169 void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects()
 170 {
 171   // Create histograms, called once.
 172   if(this->fDebug)
 173     printf("::UserCreateOutputObjects() starting\n");
 174
 175   fSelTracks = new TObjArray();
 176
 177   fSelPrimTracks = new TObjArray();
 178
 179   if (TClass::GetClass("AliPhotonHeaderObj"))
 180       TClass::GetClass("AliPhotonHeaderObj")->IgnoreTObjectStreamer();
 181   fHeader = new AliPhotonHeaderObj;
 182
 183   if (TClass::GetClass("AliPhotonConvObj"))
 184       TClass::GetClass("AliPhotonConvObj")->IgnoreTObjectStreamer();
 185   fPhotConvArray = new TClonesArray("AliPhotonConvObj");
 186
 187   if (TClass::GetClass("AliPhotonClusterObj"))
 188       TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
 189   fMyClusts = new TClonesArray("AliPhotonClusterObj");
 190
 191   if (TClass::GetClass("AliPhotonCellObj"))
 192       TClass::GetClass("AliPhotonCellObj")->IgnoreTObjectStreamer();
 193   fMyCells = new TClonesArray("AliPhotonCellObj");
 194
 195   if (TClass::GetClass("AliPhotonTrackObj"))
 196       TClass::GetClass("AliPhotonTrackObj")->IgnoreTObjectStreamer();
 197   fMyTracks = new TClonesArray("AliPhotonTrackObj");
 198
 199   if (fClusterizer || fIsTrain){
 200     if(fClusterizer)
 201       fCaloClustersName = fClusterizer->GetNewClusterArrayName();
 202     else {
 203       if(fPeriod.Contains("10h") || fPeriod.Contains("11h"))
 204         fCaloClustersName = "EmcalClustersv1";
 205       else
 206         fCaloClustersName = "EmcalClustersv2";
 207     }
 208     if (TClass::GetClass("AliPhotonClusterObj"))
 209         TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
 210     fMyAltClusts = new TClonesArray("AliPhotonClusterObj");
 211   }
 212   cout<<fCaloClustersName<<endl;
 213   if(fIsMC){
 214     if (TClass::GetClass("AliPhotonMcPartObj"))
 215         TClass::GetClass("AliPhotonMcPartObj")->IgnoreTObjectStreamer();
 216     fMyMcParts = new TClonesArray("AliPhotonMcPartObj");
 217   }
 218
 219   fCaloClusters = new TClonesArray();
 220
 221   fOutputList = new TList();
 222   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 223
 224   if( !fTree){
 225     TFile *f = OpenFile(2);
 226     TDirectory::TContext context(f);
 227     if (f) {
 228       f->SetCompressionLevel(2);
 229       fTree = new TTree("photon_ana_out", "StandaloneTree");
 230       fTree->SetDirectory(f);
 231       if (fIsTrain) {
 232         fTree->SetAutoFlush(-2*1024*1024);
 233         fTree->SetAutoSave(0);
 234       } else {
 235         fTree->SetAutoFlush(-32*1024*1024);
 236         fTree->SetAutoSave(0);
 237       }
 238
 239     fTree->Branch("header", &fHeader, 16*1024, 99);
 240     fTree->Branch("conversions", &fPhotConvArray, 8*16*1024, 99);
 241     fTree->Branch("clusters", &fMyClusts, 8*16*1024, 99);
 242     if(fClusterizer || fIsTrain)
 243       fTree->Branch(fCaloClustersName, &fMyAltClusts, 8*16*1024, 99);
 244     fTree->Branch("cells", &fMyCells, 8*16*1024, 99);
 245     fTree->Branch("IsoTracks", &fMyTracks, 8*16*1024, 99);
 246     if(fIsMC)
 247       fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99);
 248     }
 249   }
 250   //if(fIsGrid)fOutputList->Add(fTree);
 251   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 252   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
 253   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
 254
 255
 256   fNV0sBefAndAftRerun = new TH2F("hNV0sBefAndAftRerun","check if the number of v0s change with rerun;old v0 n;new v0 n",50,0.5,50.5,50,0.5,50.5);
 257   fOutputList->Add(fNV0sBefAndAftRerun);
 258
 259   fConversionVtxXY = new TH2F("hConversionVtxXY","x and y of conversion vertex candidates;x;y",1000,-100,100,1000,-100,100);
 260   fOutputList->Add(fConversionVtxXY);
 261
 262   fInvMassV0 = new TH1F("hInvMassV0","v0->GetEffMass();v0->GetEffMass();dN/dM",400,0,4);
 263   fOutputList->Add(fInvMassV0);
 264
 265   fInvMassV0KF = new TH1F("hInvMassV0KF","Inv. mass calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
 266   fOutputList->Add(fInvMassV0KF);
 267
 268   fInvMassV0SS = new TH1F("hInvMassV0SS","Inv. mass (same sign) calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
 269   fOutputList->Add(fInvMassV0SS);
 270
 271   fDedxPAll = new TH2F("hDedxPAll","dE/dx vs p (all selected tracks);p (GeV/c);dE/dx (a.u.)",400,0,40, 150,0,150);
 272   fOutputList->Add(fDedxPAll);
 273
 274
 275   PostData(1, fOutputList);
 276   PostData(2, fTree);
 277
 278   if(this->fDebug)
 279     printf("::UserCreateOutputObjects() DONE!\n");
 280
 281 }
 282
 283 //________________________________________________________________________
 284 void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *)
 285 {
 286   // User exec, called once per event.
 287
 288   Bool_t isSelected = kTRUE;
 289   if(fPeriod.Contains("11")){
 290     if(fPeriod.Contains("11a"))
 291       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 292     if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") )
 293       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 294     if(fPeriod.Contains("11h") )
 295       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);//kEMCEGA);
 296
 297   }
 298   if(fIsMC){
 299     isSelected = kTRUE;
 300     if(this->fDebug)
 301       printf("+++Message+++: MC input files.\n");
 302   }
 303   if(!isSelected){
 304     if(this->fDebug)
 305       printf("+++Message+++: Event didn't pass the selection\n");
 306     return;
 307   }
 308   if(this->fDebug){
 309     printf("::UserExec(): event accepted\n");
 310     if(fIsMC)
 311       printf("\t in MC mode\n");
 312   }
 313
 314   TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 315   TFile *inpfile = (TFile*)tree->GetCurrentFile();
 316
 317   fVev = (AliVEvent*)InputEvent();
 318   if (!fVev) {
 319     printf("ERROR: event not available\n");
 320     return;
 321   }
 322   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 323   if(this->fDebug)
 324     printf("run number = %d\n",runnumber);
 325   fESD = dynamic_cast<AliESDEvent*>(fVev);
 326   if(!fESD){
 327     fAOD = dynamic_cast<AliAODEvent*>(fVev);
 328     if(!fAOD){
 329       printf("ERROR: Invalid type of event!!!\n");
 330       return;
 331     }
 332     else if(this->fDebug)
 333       printf("AOD EVENT!\n");
 334   }
 335
 336   AliVVertex *pv = (AliVVertex*)fVev->GetPrimaryVertex();
 337   Bool_t pvStatus = kTRUE;
 338   if(fESD){
 339     AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
 340     pvStatus = esdv->GetStatus();
 341   }
 342
 343   if(!pv) {
 344     printf("Error: no primary vertex found!\n");
 345     return;
 346   }
 347   if(!pvStatus && this->fDebug)
 348     printf("bad pv status\n");
 349   if(TMath::Abs(pv->GetZ())>15)
 350     return;
 351   if(this->fDebug)
 352     printf("+++Message+++: event passed the vertex cut\n");
 353
 354   if(fESD)
 355     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 356   if(fAOD)
 357     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 358
 359   if(!fTracks){
 360     AliError("Track array in event is NULL!");
 361     if(this->fDebug)
 362       printf("returning due to not finding tracks!\n");
 363     return;
 364   }
 365
 366   const Int_t Ntracks = fTracks->GetEntriesFast();
 367   // Track loop to fill a pT spectrum
 368   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 369     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 370     if (!track)
 371       continue;
 372     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 373     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 374     if (esdTrack){
 375       if (fTrCuts && fTrCuts->IsSelected(track)){
 376         fSelTracks->Add(track);
 377         fDedxPAll->Fill(track->P(), track->GetTPCsignal());
 378       }
 379       if (fPrTrCuts && fPrTrCuts->IsSelected(track))
 380         fSelPrimTracks->Add(track);
 381     }
 382     else if(aodTrack)
 383       fSelPrimTracks->Add(track);
 384   }//track loop
 385
 386
 387
 388   fHeader->fInputFileName  = inpfile->GetName();
 389   fHeader->fRunNumber =  runnumber;
 390   fHeader->fTrClassMask    = fVev->GetHeader()->GetTriggerMask();
 391   fHeader->fTrCluster      = fVev->GetHeader()->GetTriggerCluster();
 392   AliCentrality *cent = InputEvent()->GetCentrality();
 393   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
 394   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
 395   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
 396
 397   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 398   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 399     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 400       break;
 401     /*if(fESD)
 402       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 403       else*/
 404     // if(event->GetEMCALMatrix(mod))
 405     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 406     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 407   }
 408   Int_t trackMult = 0;
 409   if(fESD){
 410     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 411     trackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 412     if(this->fDebug)
 413       printf("ESD Track mult= %d\n",trackMult);
 414   }
 415   else if(fAOD){
 416     trackMult = Ntracks;
 417     if(this->fDebug)
 418       printf("AOD Track mult= %d\n",trackMult);
 419   }
 420   if(fESD){
 421     TList *l = fESD->GetList();
 422     fCaloClusters =  dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 423     fVCells = fESD->GetEMCALCells();
 424     fHeader->fNCells = fVCells->GetNumberOfCells();
 425     if(this->fDebug)
 426       printf("ESD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
 427   }
 428   else if(fAOD){
 429     fCaloClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 430     fVCells = fAOD->GetEMCALCells();
 431     fHeader->fNCells = fVCells->GetNumberOfCells();
 432     if(this->fDebug)
 433       printf("AOD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
 434   }
 435
 436
 437   fHeader->fNClus = fCaloClusters->GetEntriesFast();
 438   fHeader->fTrackMult = trackMult;
 439
 440   fMCEvent = MCEvent();
 441   if(fMCEvent){
 442     fStack = (AliStack*)fMCEvent->Stack();
 443     if(this->fStack)
 444       fHeader->fNMcParts = this->fStack->GetNtrack();
 445     else{
 446       printf("Stack not found\n");
 447       fHeader->fNMcParts = 0;
 448       fAODMCParticles = (TClonesArray*)fVev->FindListObject("mcparticles");
 449     }
 450     if(fAODMCParticles)
 451       fHeader->fNMcParts = fAODMCParticles->GetEntriesFast();
 452   }
 453   else{
 454     if(fIsMC){
 455       printf("ERROR: MC Event not available, returning...\n");
 456       return;
 457     }
 458   }
 459
 460
 461   FindConversions();
 462   if(this->fDebug)
 463     printf("FindConversions done\n");
 464   FillMyCells();
 465   if(this->fDebug)
 466     printf("FillMyCells done\n");
 467   FillMyClusters();
 468   if(this->fDebug)
 469     printf("FillMyClusters done\n");
 470   if(fCaloClustersNew)
 471     FillMyAltClusters();
 472   FillIsoTracks();
 473   if(fIsMC)
 474     GetMcParts();
 475
 476   if(this->fDebug && fIsMC)
 477     printf("fMyMcParts nentries=%d",fMyMcParts->GetEntries());
 478
 479   fTree->Fill();
 480   fSelTracks->Clear();
 481   fSelPrimTracks->Clear();
 482   fPhotConvArray->Clear();
 483   fMyClusts->Clear();
 484   if(fMyAltClusts)
 485     fMyAltClusts->Clear();
 486   fMyCells->Clear();
 487   fMyTracks->Clear();
 488   if(fMyMcParts)
 489     fMyMcParts->Clear();
 490   fMyMcIndex = 0;
 491   fCaloClusters->Clear();
 492   if(fCaloClustersNew)
 493     fCaloClustersNew->Clear();
 494   if(fVCells)
 495     fVCells = 0;
 496   // Post output data.
 497   PostData(1, fOutputList);
 498   PostData(2, fTree);
 499 }
 500
 501 //________________________________________________________________________
 502 void AliAnalysisTaskEMCALPhoton::FindConversions() //WARNING: not ready to use with AODs
 503 {
 504   // Find conversion.
 505   if(!fESD)//not working with AODs yet
 506     return;
 507   if(this->fDebug)
 508     printf("::FindConversions() starting\n");
 509   if(!fTrCuts)
 510     return;
 511   Int_t iconv = 0;
 512   Int_t nV0Orig = 0;
 513   if(fESD)
 514     nV0Orig = fESD->GetNumberOfV0s();
 515   if(fAOD)
 516     nV0Orig = fAOD->GetNumberOfV0s();
 517   Int_t nV0s = nV0Orig;
 518   Int_t ntracks = fSelTracks->GetEntriesFast();
 519   if(fRedoV0 && !fAOD){
 520     fESD->ResetV0s();
 521     AliV0vertexer lV0vtxer;
 522     lV0vtxer.Tracks2V0vertices(fESD);
 523     nV0s = fESD->GetNumberOfV0s();
 524   }
 525   fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
 526   for(Int_t iv0=0; iv0<nV0s; iv0++){
 527     AliESDv0 * v0 = fESD->GetV0(iv0);
 528     if(!v0)
 529       continue;
 530     Double_t vpos[3];
 531     fInvMassV0->Fill(v0->GetEffMass());
 532     v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
 533     Int_t ipos = v0->GetPindex();
 534     Int_t ineg = v0->GetNindex();
 535     if(ipos<0 || ipos> ntracks)
 536       continue;
 537     if(ineg<0 || ineg> ntracks)
 538       continue;
 539     AliESDtrack *pos = static_cast<AliESDtrack*>(fSelTracks->At(ipos));
 540     AliESDtrack *neg = static_cast<AliESDtrack*>(fSelTracks->At(ineg));
 541     if(!pos || !neg)
 542       continue;
 543     /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
 544       continue;*/
 545     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
 546     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
 547     if(!paramPos || !paramNeg)
 548       continue;
 549     if(pos->GetSign() <0){//change tracks
 550       pos=neg ;
 551       neg=fESD->GetTrack(v0->GetPindex()) ;
 552       paramPos=paramNeg ;
 553       paramNeg=v0->GetParamP() ;
 554     }
 555     AliKFParticle negKF(*paramNeg,11);
 556     AliKFParticle posKF(*paramPos,-11);
 557     AliKFParticle photKF(negKF,posKF) ;
 558
 559     if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
 560       continue ;
 561
 562     if(pos->GetSign() == neg->GetSign()){
 563       fInvMassV0SS->Fill(photKF.GetMass());
 564       continue;
 565     }
 566     fInvMassV0KF->Fill(photKF.GetMass());
 567     TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
 568     if(photLV.M()>0.05 || photLV.M()<0)
 569       continue;
 570     fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
 571     Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
 572     if(convPhi<0)
 573       convPhi+=TMath::TwoPi();
 574     TVector3 vecpos(vpos);
 575     Double_t v0Phi = 0;
 576     if(vecpos.Perp()>0)
 577       vecpos.Phi();
 578     if(v0Phi<0)
 579       v0Phi+=TMath::TwoPi();
 580     AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
 581     myconv->fPt = photLV.Pt();
 582     myconv->fEta = photLV.Eta();
 583     myconv->fPhi = convPhi;
 584     myconv->fVR = vecpos.Perp();
 585     if(vecpos.Perp()>0)
 586       myconv->fVEta = vecpos.Eta();
 587     myconv->fVPhi = v0Phi;
 588     myconv->fMass = photLV.M();
 589     myconv->fMcLabel = -3; //WARNING: include the correct labeling
 590     //negative daughter
 591    myconv->fNegPt      =  negKF.GetPt();
 592    myconv->fNegEta     =  negKF.GetEta();
 593    Double_t trackPhi   =  negKF.GetPhi();
 594    if(trackPhi<0)
 595      trackPhi+=TMath::TwoPi();
 596    myconv->fNegPhi     =  trackPhi;
 597    myconv->fNegDedx    =  neg->GetTPCsignal();
 598    myconv->fNegMcLabel =  neg->GetLabel();
 599     //negative daughter
 600    myconv->fPosPt      =  posKF.GetPt();
 601    myconv->fPosEta     =  posKF.GetEta();
 602    trackPhi            =  posKF.GetPhi();
 603    if(trackPhi<0)
 604      trackPhi+=TMath::TwoPi();
 605    myconv->fPosPhi     =  trackPhi;
 606    myconv->fPosDedx    =  pos->GetTPCsignal();
 607    myconv->fPosMcLabel =  pos->GetLabel();
 608   }
 609   if(this->fDebug)
 610     printf("::FindConversions() returning...\n\n");
 611   return;
 612 }
 613
 614 //________________________________________________________________________
 615 void AliAnalysisTaskEMCALPhoton::FillMyCells()
 616 {
 617   // Fill cells.
 618   if(this->fDebug)
 619     printf("::FillMyCells() starting\n");
 620
 621   if (!fVCells)
 622     return;
 623   Int_t ncells = fVCells->GetNumberOfCells();
 624   Int_t mcel = 0;//, maxcelid=-1;
 625   Double_t maxcellE = 0;//, maxcellEta=0, maxcellPhi=0;
 626   for(Int_t icell = 0; icell<ncells; icell++){
 627     Int_t absID = TMath::Abs(fVCells->GetCellNumber(icell));
 628     AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
 629     Float_t eta=-1, phi=-1;
 630     if(!fGeom){
 631       std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
 632       return;
 633     }
 634     if(!fGeom)
 635       return;
 636     /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
 637     if(maxcellE<fVCells->GetCellAmplitude(absID)){
 638       maxcellE = fVCells->GetCellAmplitude(absID);
 639       /*maxcellEta = eta;
 640       maxcellPhi = phi;
 641       maxcelid = absID;*/
 642     }
 643     Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 644     mycell->fAbsID = absID;
 645     mycell->fE = fVCells->GetCellAmplitude(absID);
 646     mycell->fEt = fVCells->GetCellAmplitude(absID)*TMath::Sin(theta);
 647     mycell->fEta = eta;
 648     mycell->fPhi = phi;
 649     mycell->fTime = fVCells->GetCellTime(absID);
 650   }
 651   if(this->fDebug)
 652     printf("::FillMyCells() returning...\n\n");
 653 }
 654
 655 //________________________________________________________________________
 656 void AliAnalysisTaskEMCALPhoton::FillMyClusters()
 657 {
 658   // Fill clusters.
 659   if(this->fDebug)
 660     printf("::FillMyClusters() starting\n");
 661
 662   if (!fCaloClusters){
 663     printf("CaloClusters is empty!\n");
 664     return;
 665   }
 666   Int_t nclus = fCaloClusters->GetEntries();
 667   if(0==nclus)
 668     printf("CaloClusters has ZERO entries\n");
 669   Int_t mcl = 0, maxcelid=-1;
 670   Double_t maxcellE=0, maxcellEtac=0,maxcellPhic=0;
 671   for(Int_t ic=0; ic < nclus; ic++){
 672     AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(ic));
 673     if(!clus)
 674       continue;
 675     if(!clus->IsEMCAL())
 676       continue;
 677     if(clus->E() < fClusThresh)
 678       continue;
 679     if(fDebug)
 680       printf("cluster %d survived\n", ic);
 681     Float_t pos[3];
 682     clus->GetPosition(pos);
 683     TVector3 cpos(pos);
 684     TString cellsAbsID;
 685     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
 686     myclus->fE       = clus->E();
 687     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 688     myclus->fR       = cpos.Perp();
 689     myclus->fEta     = cpos.Eta();
 690     myclus->fPhi     = cpos.Phi();
 691     if(cpos.Phi()<0){
 692       myclus->fPhi+=TMath::TwoPi();
 693     }
 694     myclus->fN       = clus->GetNCells();
 695     Short_t  id = -1;
 696     myclus->fEmax    = GetMaxCellEnergy( clus, id);
 697     myclus->fIdmax   = id;
 698     if(maxcellE <  myclus->fEmax){
 699       maxcellE =  myclus->fEmax;
 700       maxcelid = id;
 701       maxcellEtac = cpos.Eta();
 702       maxcellPhic = cpos.Phi();
 703     }
 704     myclus->fTmax    = fVCells->GetCellTime(id);
 705     myclus->fEcross  = GetCrossEnergy( clus, id);
 706     myclus->fDisp    = clus->GetDispersion();
 707     myclus->fM20     = clus->GetM20();
 708     myclus->fM02     = clus->GetM02();
 709     myclus->fTrDEta  = clus->GetTrackDz();
 710     myclus->fTrDPhi  = clus->GetTrackDx();
 711     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 712     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 713     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 714     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 715     myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 716     myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 717     myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 718     myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 719     for(Int_t icell=0;icell<myclus->fN;icell++){
 720       Int_t absID = clus->GetCellAbsId(icell);
 721       cellsAbsID.Append(Form("%d",absID));
 722       cellsAbsID.Append(";");
 723     }
 724     myclus->fCellsAbsId = cellsAbsID;
 725     myclus->fMcLabel = clus->GetLabel();
 726     Int_t matchIndex = clus->GetTrackMatchedIndex();
 727     if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
 728       myclus->fTrEp = -1;
 729       continue;
 730     }
 731     AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
 732     if(!track)
 733       continue;
 734     if(fESD){
 735       if(!fPrTrCuts)
 736         continue;
 737       if(!fPrTrCuts->IsSelected(track))
 738         continue;
 739     }
 740
 741     myclus->fTrEp = clus->E()/track->P();
 742     myclus->fTrDedx = track->GetTPCsignal();
 743   }
 744   if(this->fDebug){
 745     printf(" ---===+++ Max Cell among clusters: id=%d, E=%1.2f, eta-clus=%1.2f, phi-clus=%1.2f\n",maxcelid,maxcellE,maxcellEtac,maxcellPhic);
 746     printf("::FillMyClusters() returning...\n\n");
 747   }
 748
 749 }
 750 //________________________________________________________________________
 751 void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
 752 {
 753   // Fill clusters.
 754   if(this->fDebug)
 755     printf("::FillMyAltClusters() starting\n");
 756
 757   if(!fCaloClustersNew)
 758     return;
 759   Int_t nclus = fCaloClustersNew->GetEntries();
 760   Int_t mcl = 0;
 761   for(Int_t ic=0; ic < nclus; ic++){
 762     AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic));
 763     if(!clus)
 764       continue;
 765     if(!clus->IsEMCAL())
 766       continue;
 767     if(clus->E() < fClusThresh)
 768       continue;
 769     Float_t pos[3];
 770     clus->GetPosition(pos);
 771     TVector3 cpos(pos);
 772     TString cellsAbsID;
 773     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
 774     myclus->fE       = clus->E();
 775     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 776     myclus->fR       = cpos.Perp();
 777     myclus->fEta     = cpos.Eta();
 778     myclus->fPhi     = cpos.Phi();
 779     if(cpos.Phi()<0){
 780       myclus->fPhi+=TMath::TwoPi();
 781     }
 782     myclus->fN       = clus->GetNCells();
 783     myclus->fDisp    = clus->GetDispersion();
 784     myclus->fM20     = clus->GetM20();
 785     myclus->fM02     = clus->GetM02();
 786     myclus->fTrDEta  = clus->GetTrackDz();
 787     myclus->fTrDPhi  = clus->GetTrackDx();
 788     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 789     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 790     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 791     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 792     for(Int_t icell=0;icell<myclus->fN;icell++){
 793       Int_t absID = clus->GetCellAbsId(icell);
 794       cellsAbsID.Append(Form("%d",absID));
 795       cellsAbsID.Append(";");
 796     }
 797     myclus->fCellsAbsId = cellsAbsID;
 798     myclus->fMcLabel = clus->GetLabel();
 799     Int_t matchIndex = clus->GetTrackMatchedIndex();
 800     if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
 801       myclus->fTrEp = -1;
 802       continue;
 803     }
 804     AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
 805     if(!track){
 806       myclus->fTrEp = -1;
 807       continue;
 808     }
 809     if(!fPrTrCuts){
 810       myclus->fTrEp = -1;
 811       continue;
 812     }
 813     if(!fPrTrCuts->IsSelected(track)){
 814       myclus->fTrEp = -1;
 815       continue;
 816     }
 817     myclus->fTrEp = clus->E()/track->P();
 818   }
 819   if(this->fDebug)
 820     printf("::FillMyAltClusters() returning...\n\n");
 821
 822 }
 823 //________________________________________________________________________
 824 void  AliAnalysisTaskEMCALPhoton::FillIsoTracks()
 825 {
 826   // Fill high pt tracks.
 827   if(this->fDebug)
 828     printf("::FillIsoTracks() starting\n");
 829
 830   if(!fSelPrimTracks)
 831     return;
 832   Int_t ntracks = fSelPrimTracks->GetEntries();
 833   Int_t imtr = 0;
 834   for(Int_t it=0;it<ntracks; it++){
 835     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it));
 836     if(!track)
 837       continue;
 838     /*if(track->Phi()<1.0 || track->Phi()>3.55)
 839       continue;*/
 840     if(TMath::Abs(track->Eta())>1.1)
 841       continue;
 842     /*if(track->Pt()<3)
 843       continue;*/
 844     AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
 845     mtr->fPt = track->Pt();
 846     mtr->fEta = track->Eta();
 847     mtr->fPhi = track->Phi();
 848     mtr->fCharge = track->Charge();
 849     mtr->fDedx = track->GetTPCsignal();
 850     mtr->fMcLabel = track->GetLabel();
 851   }
 852   if(this->fDebug)
 853     printf("::FillIsoTracks() returning...\n\n");
 854 }
 855
 856 //________________________________________________________________________
 857 void  AliAnalysisTaskEMCALPhoton::GetMcParts()
 858 {
 859    // Get MC particles.
 860   if(this->fDebug)
 861     printf("::GetMcParts() starting\n");
 862
 863   if (!this->fStack && !fAODMCParticles)
 864     return;
 865   if(this->fDebug)
 866     printf("either stack or aodmcpaticles exists\n");
 867   const AliVVertex *evtVtx = 0;
 868   if(this->fStack)
 869      evtVtx = fMCEvent->GetPrimaryVertex();
 870   else
 871     printf("no such thing as mc vvtx\n");
 872   /*if (!evtVtx)
 873     return;*/
 874   if(this->fDebug)
 875     printf("mc vvertex ok\n");
 876   Int_t nTracks = 0;
 877   if(this->fStack)
 878     nTracks = this->fStack->GetNtrack();
 879   else if(fAODMCParticles)
 880     nTracks = fAODMCParticles->GetEntriesFast();
 881   TParticle *mcP = 0;
 882   AliAODMCParticle *amcP = 0;
 883   if(this->fDebug)
 884     printf("loop in the %d mc particles starting\n",nTracks);
 885   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
 886     if(this->fStack)
 887       mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack));
 888     if(fAODMCParticles)
 889       amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
 890
 891     // primary particle
 892     Double_t dR = 0;
 893     if(mcP){
 894       dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
 895                        (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
 896                        (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
 897     }
 898     if((dR > 0.5))
 899       continue;
 900
 901
 902     // kinematic cuts
 903     Double_t pt = 0;
 904     Double_t eta = 0;
 905     Double_t phi  = 0;
 906     Int_t mother = -1;
 907     Int_t pdgcode = 0;
 908      if(mcP){
 909       pt = mcP->Pt() ;
 910       eta = mcP->Eta();
 911       phi = mcP->Phi();
 912       mother = mcP->GetMother(0);
 913       pdgcode = mcP->GetPdgCode();
 914      } else {
 915        if(amcP){
 916          pt = amcP->Pt();
 917          eta = amcP->Eta();
 918          phi = amcP->Phi();
 919          mother = amcP->GetMother();
 920          pdgcode = amcP->GetPdgCode();
 921        } else
 922        continue;
 923      }
 924     if (pt<0.5)
 925       continue;
 926
 927     if (TMath::Abs(eta)>0.7)
 928       continue;
 929
 930     if (phi<1.0||phi>3.3)
 931       continue;
 932
 933     if (mother!=6 && mother!=7 )
 934       continue;
 935
 936
 937     // pion or eta meson or direct photon
 938     if(pdgcode == 111) {
 939     } else if(pdgcode == 221) {
 940     } else if(pdgcode == 22 ) {
 941     } else {
 942       continue;
 943     }
 944
 945     FillMcPart( fMyMcIndex++, iTrack);
 946   }
 947   if(this->fDebug)
 948     printf("::GetMcParts() returning...\n\n");
 949 }
 950
 951 //________________________________________________________________________
 952 void  AliAnalysisTaskEMCALPhoton::FillMcPart(  Int_t itrack, Int_t label)
 953 {
 954   // Fill MC particles.
 955   Int_t nTracks = 0;
 956   TParticle *mcP = 0;
 957   AliAODMCParticle *amcP= 0;
 958   if(this->fStack){
 959     nTracks = this->fStack->GetNtrack();
 960     mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
 961   }
 962   else if(fAODMCParticles){
 963     nTracks = fAODMCParticles->GetEntriesFast();
 964     amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack));
 965  }
 966   if(this->fDebug)
 967     printf("\t::FillMcParts() starting with label %d\n", itrack);
 968    TVector3 vmcv;
 969    if(mcP)
 970      vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz());
 971    else {
 972      if(amcP)
 973        vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv());
 974      else
 975        return;
 976    }
 977
 978   AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
 979   mcp->fLabel = label ;
 980   if(mcP){
 981     mcp->fPdg = mcP->GetPdgCode() ;
 982     mcp->fPt = mcP->Pt() ;
 983     mcp->fEta = mcP->Eta() ;
 984     mcp->fPhi = mcP->Phi() ;
 985     mcp->fMother = mcP->GetMother(0) ;
 986     mcp->fFirstD = mcP->GetFirstDaughter() ;
 987     mcp->fLastD = mcP->GetLastDaughter() ;
 988     mcp->fStatus = mcP->GetStatusCode();
 989   }
 990   if(amcP){
 991     mcp->fPdg = amcP->GetPdgCode() ;
 992     mcp->fPt = amcP->Pt() ;
 993     mcp->fEta = amcP->Eta() ;
 994     mcp->fPhi = amcP->Phi() ;
 995     mcp->fMother = amcP->GetMother() ;
 996     mcp->fFirstD = amcP->GetDaughter(0) ;
 997     mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ;
 998     mcp->fStatus = amcP->GetStatus();
 999   }
1000   mcp->fVR = vmcv.Perp();
1001   if(vmcv.Perp()>0){
1002     mcp->fVEta = vmcv.Eta() ;
1003     mcp->fVPhi = vmcv.Phi() ;
1004   }
1005   if(itrack == 8){
1006     mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2);
1007     mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2);
1008     }
1009   if(this->fDebug)
1010     printf("\t::FillMcParts(): label=%d, pdg=%d, pt=%1.1f, mom=%d, 1stD=%d,last=%d\n\t::FillMcParts() returning...\n\n", mcp->fLabel,mcp->fPdg,mcp->fPt,mcp->fMother,mcp->fFirstD,mcp->fLastD);
1011   for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){
1012     if(id<=mcp->fMother)
1013       continue;
1014     if(id<0 || id>nTracks)
1015       continue;
1016     FillMcPart( fMyMcIndex++, id);
1017   }
1018
1019 }
1020 //________________________________________________________________________
1021 Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const
1022 {
1023   if(this->fDebug){
1024     printf("\t\t::GetMcIsolation() starting\n");
1025     //printf("\t\t   incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack);
1026   }
1027   if (!this->fStack && !this->fAODMCParticles && this->fDebug){
1028     printf("\t\t\tNo MC stack/array!\n");
1029     return -1;
1030   }
1031   if(itrack<6 || itrack>8)
1032     return -1;
1033   if(this->fDebug)
1034     printf("\t\t\tparticle of interest selected\n");
1035   TParticle *mcP = 0;
1036   AliAODMCParticle *amcP = 0;
1037   Int_t pdgcode = 0;
1038   Float_t eta = 0;
1039   Float_t phi = 0;
1040   Double_t sumpt=0;
1041   Float_t dR;
1042   Int_t nparts =  0;
1043   if(this->fStack){
1044     nparts = fStack->GetNtrack();
1045     mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
1046     eta = mcP->Eta();
1047     phi = mcP->Phi();
1048     pdgcode = mcP->GetPdgCode();
1049   }
1050   if(this->fAODMCParticles){
1051     nparts = fAODMCParticles->GetEntriesFast();
1052     amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack));
1053     if(amcP){
1054       eta = amcP->Eta();
1055       phi = amcP->Phi();
1056       pdgcode = amcP->GetPdgCode();
1057     }
1058   }
1059   if(pdgcode!=22)
1060     return -1;
1061   TParticle *mcisop = 0;
1062   AliAODMCParticle *amcisop = 0;
1063   for(Int_t ip = 0; ip<nparts; ip++){
1064     if(ip==itrack)
1065       continue;
1066     if(this->fStack)
1067       mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip));
1068     if(fAODMCParticles)
1069       amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip));
1070     Int_t status = 0;
1071     Int_t mother = 0;
1072     Float_t pt = 0;
1073     Float_t isophi = -99;
1074     Float_t isoeta = -99;
1075     TVector3 vmcv;
1076     if(mcisop){
1077       status = mcisop->GetStatusCode();
1078       mother = mcisop->GetMother(0);
1079       pt = mcisop->Pt();
1080       isophi = mcisop->Phi();
1081       isoeta = mcisop->Eta();
1082       vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
1083     }
1084     else {
1085       if(amcisop){
1086         status = amcisop->GetStatus();
1087         mother = amcisop->GetMother();
1088         pt = amcisop->Pt();
1089         isophi = amcisop->Phi();
1090         isoeta = amcisop->Eta();
1091         vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv());
1092       }
1093       else
1094         continue;
1095     }
1096     if(status!=1)
1097       continue;
1098     if(mother == itrack)
1099       continue;
1100     if(pt<ptcut)
1101       continue;
1102     if(vmcv.Perp()>1)
1103       continue;
1104     dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta));
1105     if(dR>radius)
1106       continue;
1107     sumpt += pt;
1108   }
1109   if(this->fDebug)
1110     printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt);
1111   return sumpt;
1112  }
1113
1114 //________________________________________________________________________
1115 Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1116 {
1117   // Compute isolation based on tracks.
1118   if(this->fDebug)
1119     printf("\t::GetTrackIsolation() starting\n");
1120
1121   Double_t trkIsolation = 0;
1122   Double_t rad2 = radius*radius;
1123   Int_t ntrks = fSelPrimTracks->GetEntries();
1124   for(Int_t j = 0; j<ntrks; ++j) {
1125     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
1126     if (!track)
1127       continue;
1128     if (track->Pt()<pt)
1129       continue;
1130
1131     Float_t eta = track->Eta();
1132     Float_t phi = track->Phi();
1133     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1134     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1135     if(dist>rad2)
1136       continue;
1137     trkIsolation += track->Pt();
1138   }
1139   if(this->fDebug)
1140     printf("\t::GetTrackIsolation() returning\n\n");
1141   return trkIsolation;
1142 }
1143
1144 //________________________________________________________________________
1145 Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
1146 {
1147   // Get phi band.
1148
1149   if(!fSelPrimTracks)
1150     return 0;
1151   Int_t ntracks = fSelPrimTracks->GetEntries();
1152   Double_t loweta = eta - R;
1153   Double_t upeta = eta + R;
1154   Double_t upphi = phi + R;
1155   Double_t lowphi = phi - R;
1156   Double_t et = 0;
1157   for(Int_t itr=0; itr<ntracks; itr++){
1158     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr));
1159     if(!track)
1160       continue;
1161     if(track->Pt()<minpt)
1162       continue;
1163     if((track->Eta() < upeta) && (track->Eta() > loweta))
1164       continue;
1165     if((track->Phi() > upphi) || (track->Phi() < lowphi))
1166       continue;
1167     et+=track->Pt();
1168   }
1169   return et;
1170 }
1171
1172 //________________________________________________________________________
1173 Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1174 {
1175   // Calculate the energy of cross cells around the leading cell.
1176   if(!fVCells)
1177     return 0;
1178
1179   if (!fGeom)
1180     return 0;
1181
1182   Int_t iSupMod = -1;
1183   Int_t iTower  = -1;
1184   Int_t iIphi   = -1;
1185   Int_t iIeta   = -1;
1186   Int_t iphi    = -1;
1187   Int_t ieta    = -1;
1188   Int_t iphis   = -1;
1189   Int_t ietas   = -1;
1190
1191   Double_t crossEnergy = 0;
1192
1193   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1194   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1195
1196   Int_t ncells = cluster->GetNCells();
1197   for (Int_t i=0; i<ncells; i++) {
1198     Int_t cellAbsId = cluster->GetCellAbsId(i);
1199     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1200     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1201     Int_t aphidiff = TMath::Abs(iphi-iphis);
1202     if (aphidiff>1)
1203       continue;
1204     Int_t aetadiff = TMath::Abs(ieta-ietas);
1205     if (aetadiff>1)
1206       continue;
1207     if ( (aphidiff==1 && aetadiff==0) ||
1208         (aphidiff==0 && aetadiff==1) ) {
1209       crossEnergy += fVCells->GetCellAmplitude(cellAbsId);
1210     }
1211   }
1212
1213   return crossEnergy;
1214 }
1215
1216 //________________________________________________________________________
1217 Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1218 {
1219   // Get maximum energy of attached cell.
1220
1221   id = -1;
1222   if(!fVCells)
1223     return 0;
1224
1225   Double_t maxe = 0;
1226   Int_t ncells = cluster->GetNCells();
1227   for (Int_t i=0; i<ncells; i++) {
1228     Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1229     if (e>maxe) {
1230       maxe = e;
1231       id   = cluster->GetCellAbsId(i);
1232     }
1233   }
1234   return maxe;
1235 }
1236
1237 //________________________________________________________________________
1238 void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
1239 {
1240   // Called once at the end of the query
1241 /*  if(fIsGrid)
1242     return;*/
1243   if (fTree) {
1244     printf("***tree %s being saved***\n",fTree->GetName());
1245     TFile *f = OpenFile(2);
1246     TDirectory::TContext context(f);
1247     if (f)
1248       fTree->Write();
1249   }
1250 }