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(TMath::Abs(pv->GetZ())>15)
 348     return;
 349   if(this->fDebug)
 350     printf("+++Message+++: event passed the vertex cut\n");
 351
 352   if(fESD)
 353     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 354   if(fAOD)
 355     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 356
 357   if(!fTracks){
 358     AliError("Track array in event is NULL!");
 359     if(this->fDebug)
 360       printf("returning due to not finding tracks!\n");
 361     return;
 362   }
 363
 364   const Int_t Ntracks = fTracks->GetEntriesFast();
 365   // Track loop to fill a pT spectrum
 366   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 367     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 368     if (!track)
 369       continue;
 370     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 371     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 372     if (esdTrack){
 373       if (fTrCuts && fTrCuts->IsSelected(track)){
 374         fSelTracks->Add(track);
 375         fDedxPAll->Fill(track->P(), track->GetTPCsignal());
 376       }
 377       if (fPrTrCuts && fPrTrCuts->IsSelected(track))
 378         fSelPrimTracks->Add(track);
 379     }
 380     else if(aodTrack)
 381       fSelPrimTracks->Add(track);
 382   }//track loop
 383
 384
 385
 386   fHeader->fInputFileName  = inpfile->GetName();
 387   fHeader->fTrClassMask    = fVev->GetHeader()->GetTriggerMask();
 388   fHeader->fTrCluster      = fVev->GetHeader()->GetTriggerCluster();
 389   AliCentrality *cent = InputEvent()->GetCentrality();
 390   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
 391   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
 392   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
 393
 394   if(fDebug)
 395     printf("AliEMCALgeo file found\n");
 396   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 397   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 398     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 399       break;
 400     /*if(fESD)
 401       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 402       else*/
 403     // if(event->GetEMCALMatrix(mod))
 404     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 405     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 406   }
 407   Int_t trackMult = 0;
 408   if(fESD){
 409     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 410     trackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 411     if(this->fDebug)
 412       printf("ESD Track mult= %d\n",trackMult);
 413   }
 414   else if(fAOD){
 415     trackMult = Ntracks;
 416     if(this->fDebug)
 417       printf("AOD Track mult= %d\n",trackMult);
 418   }
 419   if(fESD){
 420     TList *l = fESD->GetList();
 421     fCaloClusters =  dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 422     fVCells = fESD->GetEMCALCells();
 423     fHeader->fNCells = fVCells->GetNumberOfCells();
 424     if(this->fDebug)
 425       printf("ESD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
 426   }
 427   else if(fAOD){
 428     fCaloClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 429     fVCells = fAOD->GetEMCALCells();
 430     fHeader->fNCells = fVCells->GetNumberOfCells();
 431     if(this->fDebug)
 432       printf("AOD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
 433   }
 434
 435
 436   fHeader->fNClus = fCaloClusters->GetEntriesFast();
 437   fHeader->fTrackMult = trackMult;
 438
 439   fMCEvent = MCEvent();
 440   if(fMCEvent){
 441     fStack = (AliStack*)fMCEvent->Stack();
 442     if(this->fStack)
 443       fHeader->fNMcParts = this->fStack->GetNtrack();
 444     else{
 445       printf("Stack not found\n");
 446       fHeader->fNMcParts = 0;
 447       fAODMCParticles = (TClonesArray*)fVev->FindListObject("mcparticles");
 448     }
 449     if(fAODMCParticles)
 450       fHeader->fNMcParts = fAODMCParticles->GetEntriesFast();
 451   }
 452   else{
 453     if(fIsMC){
 454       printf("ERROR: MC Event not available, returning...\n");
 455       return;
 456     }
 457   }
 458
 459
 460   FindConversions();
 461   if(this->fDebug)
 462     printf("FindConversions done\n");
 463   FillMyCells();
 464   if(this->fDebug)
 465     printf("FillMyCells done\n");
 466   FillMyClusters();
 467   if(this->fDebug)
 468     printf("FillMyClusters done\n");
 469   if(fCaloClustersNew)
 470     FillMyAltClusters();
 471   FillIsoTracks();
 472   if(fIsMC)
 473     GetMcParts();
 474
 475   if(this->fDebug && fIsMC)
 476     printf("fMyMcParts nentries=%d",fMyMcParts->GetEntries());
 477
 478   fTree->Fill();
 479   fSelTracks->Clear();
 480   fSelPrimTracks->Clear();
 481   fPhotConvArray->Clear();
 482   fMyClusts->Clear();
 483   if(fMyAltClusts)
 484     fMyAltClusts->Clear();
 485   fMyCells->Clear();
 486   fMyTracks->Clear();
 487   if(fMyMcParts)
 488     fMyMcParts->Clear();
 489   fMyMcIndex = 0;
 490   fCaloClusters->Clear();
 491   if(fCaloClustersNew)
 492     fCaloClustersNew->Clear();
 493   if(fVCells)
 494     fVCells = 0;
 495   // Post output data.
 496   PostData(1, fOutputList);
 497   PostData(2, fTree);
 498 }
 499
 500 //________________________________________________________________________
 501 void AliAnalysisTaskEMCALPhoton::FindConversions() //WARNING: not ready to use with AODs
 502 {
 503   // Find conversion.
 504   if(!fESD)//not working with AODs yet
 505     return;
 506   if(this->fDebug)
 507     printf("::FindConversions() starting\n");
 508   if(!fTrCuts)
 509     return;
 510   Int_t iconv = 0;
 511   Int_t nV0Orig = 0;
 512   if(fESD)
 513     nV0Orig = fESD->GetNumberOfV0s();
 514   if(fAOD)
 515     nV0Orig = fAOD->GetNumberOfV0s();
 516   Int_t nV0s = nV0Orig;
 517   Int_t ntracks = fSelTracks->GetEntriesFast();
 518   if(fRedoV0 && !fAOD){
 519     fESD->ResetV0s();
 520     AliV0vertexer lV0vtxer;
 521     lV0vtxer.Tracks2V0vertices(fESD);
 522     nV0s = fESD->GetNumberOfV0s();
 523   }
 524   fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
 525   for(Int_t iv0=0; iv0<nV0s; iv0++){
 526     AliESDv0 * v0 = fESD->GetV0(iv0);
 527     if(!v0)
 528       continue;
 529     Double_t vpos[3];
 530     fInvMassV0->Fill(v0->GetEffMass());
 531     v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
 532     Int_t ipos = v0->GetPindex();
 533     Int_t ineg = v0->GetNindex();
 534     if(ipos<0 || ipos> ntracks)
 535       continue;
 536     if(ineg<0 || ineg> ntracks)
 537       continue;
 538     AliESDtrack *pos = static_cast<AliESDtrack*>(fSelTracks->At(ipos));
 539     AliESDtrack *neg = static_cast<AliESDtrack*>(fSelTracks->At(ineg));
 540     if(!pos || !neg)
 541       continue;
 542     /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
 543       continue;*/
 544     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
 545     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
 546     if(!paramPos || !paramNeg)
 547       continue;
 548     if(pos->GetSign() <0){//change tracks
 549       pos=neg ;
 550       neg=fESD->GetTrack(v0->GetPindex()) ;
 551       paramPos=paramNeg ;
 552       paramNeg=v0->GetParamP() ;
 553     }
 554     AliKFParticle negKF(*paramNeg,11);
 555     AliKFParticle posKF(*paramPos,-11);
 556     AliKFParticle photKF(negKF,posKF) ;
 557
 558     if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
 559       continue ;
 560
 561     if(pos->GetSign() == neg->GetSign()){
 562       fInvMassV0SS->Fill(photKF.GetMass());
 563       continue;
 564     }
 565     fInvMassV0KF->Fill(photKF.GetMass());
 566     TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
 567     if(photLV.M()>0.05 || photLV.M()<0)
 568       continue;
 569     fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
 570     Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
 571     if(convPhi<0)
 572       convPhi+=TMath::TwoPi();
 573     TVector3 vecpos(vpos);
 574     Double_t v0Phi = 0;
 575     if(vecpos.Perp()>0)
 576       vecpos.Phi();
 577     if(v0Phi<0)
 578       v0Phi+=TMath::TwoPi();
 579     AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
 580     myconv->fPt = photLV.Pt();
 581     myconv->fEta = photLV.Eta();
 582     myconv->fPhi = convPhi;
 583     myconv->fVR = vecpos.Perp();
 584     if(vecpos.Perp()>0)
 585       myconv->fVEta = vecpos.Eta();
 586     myconv->fVPhi = v0Phi;
 587     myconv->fMass = photLV.M();
 588     myconv->fMcLabel = -3; //WARNING: include the correct labeling
 589     //negative daughter
 590    myconv->fNegPt      =  negKF.GetPt();
 591    myconv->fNegEta     =  negKF.GetEta();
 592    Double_t trackPhi   =  negKF.GetPhi();
 593    if(trackPhi<0)
 594      trackPhi+=TMath::TwoPi();
 595    myconv->fNegPhi     =  trackPhi;
 596    myconv->fNegDedx    =  neg->GetTPCsignal();
 597    myconv->fNegMcLabel =  neg->GetLabel();
 598     //negative daughter
 599    myconv->fPosPt      =  posKF.GetPt();
 600    myconv->fPosEta     =  posKF.GetEta();
 601    trackPhi            =  posKF.GetPhi();
 602    if(trackPhi<0)
 603      trackPhi+=TMath::TwoPi();
 604    myconv->fPosPhi     =  trackPhi;
 605    myconv->fPosDedx    =  pos->GetTPCsignal();
 606    myconv->fPosMcLabel =  pos->GetLabel();
 607   }
 608   if(this->fDebug)
 609     printf("::FindConversions() returning...\n\n");
 610   return;
 611 }
 612
 613 //________________________________________________________________________
 614 void AliAnalysisTaskEMCALPhoton::FillMyCells()
 615 {
 616   // Fill cells.
 617   if(this->fDebug)
 618     printf("::FillMyCells() starting\n");
 619
 620   if (!fVCells)
 621     return;
 622   Int_t ncells = fVCells->GetNumberOfCells();
 623   Int_t mcel = 0;
 624   for(Int_t icell = 0; icell<ncells; icell++){
 625     Int_t absID = TMath::Abs(fVCells->GetCellNumber(icell));
 626     AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
 627     Float_t eta=-1, phi=-1;
 628     if(!fGeom){
 629       std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
 630       return;
 631     }
 632     if(!fGeom)
 633       return;
 634     /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
 635     Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 636     mycell->fAbsID = absID;
 637     mycell->fE = fVCells->GetCellAmplitude(absID);
 638     mycell->fEt = fVCells->GetCellAmplitude(absID)*TMath::Sin(theta);
 639     mycell->fEta = eta;
 640     mycell->fPhi = phi;
 641     mycell->fTime = fVCells->GetCellTime(absID);
 642   }
 643   if(this->fDebug)
 644     printf("::FillMyCells() returning...\n\n");
 645 }
 646
 647 //________________________________________________________________________
 648 void AliAnalysisTaskEMCALPhoton::FillMyClusters()
 649 {
 650   // Fill clusters.
 651   if(this->fDebug)
 652     printf("::FillMyClusters() starting\n");
 653
 654   if (!fCaloClusters){
 655     printf("CaloClusters is empty!\n");
 656     return;
 657   }
 658   Int_t nclus = fCaloClusters->GetEntries();
 659   if(0==nclus)
 660     printf("CaloClusters has ZERO entries\n");
 661   Int_t mcl = 0;
 662   for(Int_t ic=0; ic < nclus; ic++){
 663     AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(ic));
 664     if(!clus)
 665       continue;
 666     if(!clus->IsEMCAL())
 667       continue;
 668     if(clus->E() < fClusThresh)
 669       continue;
 670     Float_t pos[3];
 671     clus->GetPosition(pos);
 672     TVector3 cpos(pos);
 673     TString cellsAbsID;
 674     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
 675     myclus->fE       = clus->E();
 676     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 677     myclus->fR       = cpos.Perp();
 678     myclus->fEta     = cpos.Eta();
 679     myclus->fPhi     = cpos.Phi();
 680     if(cpos.Phi()<0){
 681       myclus->fPhi+=TMath::TwoPi();
 682     }
 683     myclus->fN       = clus->GetNCells();
 684     Short_t  id = -1;
 685     myclus->fEmax    = GetMaxCellEnergy( clus, id);
 686     myclus->fIdmax   = id;
 687     myclus->fTmax    = fVCells->GetCellTime(id);
 688     myclus->fEcross  = GetCrossEnergy( clus, id);
 689     myclus->fDisp    = clus->GetDispersion();
 690     myclus->fM20     = clus->GetM20();
 691     myclus->fM02     = clus->GetM02();
 692     myclus->fTrDEta  = clus->GetTrackDz();
 693     myclus->fTrDPhi  = clus->GetTrackDx();
 694     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 695     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 696     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 697     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 698     myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 699     myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 700     myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 701     myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 702     for(Int_t icell=0;icell<myclus->fN;icell++){
 703       Int_t absID = clus->GetCellAbsId(icell);
 704       cellsAbsID.Append(Form("%d",absID));
 705       cellsAbsID.Append(";");
 706     }
 707     myclus->fCellsAbsId = cellsAbsID;
 708     myclus->fMcLabel = clus->GetLabel();
 709     Int_t matchIndex = clus->GetTrackMatchedIndex();
 710     if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
 711       myclus->fTrEp = -1;
 712       continue;
 713     }
 714     AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
 715     if(!track)
 716       continue;
 717     if(fESD){
 718       if(!fPrTrCuts)
 719         continue;
 720       if(!fPrTrCuts->IsSelected(track))
 721         continue;
 722     }
 723
 724     myclus->fTrEp = clus->E()/track->P();
 725     myclus->fTrDedx = track->GetTPCsignal();
 726   }
 727   if(this->fDebug)
 728     printf("::FillMyClusters() returning...\n\n");
 729
 730 }
 731 //________________________________________________________________________
 732 void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
 733 {
 734   // Fill clusters.
 735   if(this->fDebug)
 736     printf("::FillMyAltClusters() starting\n");
 737
 738   if(!fCaloClustersNew)
 739     return;
 740   Int_t nclus = fCaloClustersNew->GetEntries();
 741   Int_t mcl = 0;
 742   for(Int_t ic=0; ic < nclus; ic++){
 743     AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic));
 744     if(!clus)
 745       continue;
 746     if(!clus->IsEMCAL())
 747       continue;
 748     if(clus->E() < fClusThresh)
 749       continue;
 750     Float_t pos[3];
 751     clus->GetPosition(pos);
 752     TVector3 cpos(pos);
 753     TString cellsAbsID;
 754     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
 755     myclus->fE       = clus->E();
 756     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 757     myclus->fR       = cpos.Perp();
 758     myclus->fEta     = cpos.Eta();
 759     myclus->fPhi     = cpos.Phi();
 760     if(cpos.Phi()<0){
 761       myclus->fPhi+=TMath::TwoPi();
 762     }
 763     myclus->fN       = clus->GetNCells();
 764     myclus->fDisp    = clus->GetDispersion();
 765     myclus->fM20     = clus->GetM20();
 766     myclus->fM02     = clus->GetM02();
 767     myclus->fTrDEta  = clus->GetTrackDz();
 768     myclus->fTrDPhi  = clus->GetTrackDx();
 769     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 770     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 771     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 772     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 773     for(Int_t icell=0;icell<myclus->fN;icell++){
 774       Int_t absID = clus->GetCellAbsId(icell);
 775       cellsAbsID.Append(Form("%d",absID));
 776       cellsAbsID.Append(";");
 777     }
 778     myclus->fCellsAbsId = cellsAbsID;
 779     myclus->fMcLabel = clus->GetLabel();
 780     Int_t matchIndex = clus->GetTrackMatchedIndex();
 781     if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
 782       myclus->fTrEp = -1;
 783       continue;
 784     }
 785     AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
 786     if(!track){
 787       myclus->fTrEp = -1;
 788       continue;
 789     }
 790     if(!fPrTrCuts){
 791       myclus->fTrEp = -1;
 792       continue;
 793     }
 794     if(!fPrTrCuts->IsSelected(track)){
 795       myclus->fTrEp = -1;
 796       continue;
 797     }
 798     myclus->fTrEp = clus->E()/track->P();
 799   }
 800   if(this->fDebug)
 801     printf("::FillMyAltClusters() returning...\n\n");
 802
 803 }
 804 //________________________________________________________________________
 805 void  AliAnalysisTaskEMCALPhoton::FillIsoTracks()
 806 {
 807   // Fill high pt tracks.
 808   if(this->fDebug)
 809     printf("::FillIsoTracks() starting\n");
 810
 811   if(!fSelPrimTracks)
 812     return;
 813   Int_t ntracks = fSelPrimTracks->GetEntries();
 814   Int_t imtr = 0;
 815   for(Int_t it=0;it<ntracks; it++){
 816     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it));
 817     if(!track)
 818       continue;
 819     /*if(track->Phi()<1.0 || track->Phi()>3.55)
 820       continue;*/
 821     if(TMath::Abs(track->Eta())>1.1)
 822       continue;
 823     /*if(track->Pt()<3)
 824       continue;*/
 825     AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
 826     mtr->fPt = track->Pt();
 827     mtr->fEta = track->Eta();
 828     mtr->fPhi = track->Phi();
 829     mtr->fCharge = track->Charge();
 830     mtr->fDedx = track->GetTPCsignal();
 831     mtr->fMcLabel = track->GetLabel();
 832   }
 833   if(this->fDebug)
 834     printf("::FillIsoTracks() returning...\n\n");
 835 }
 836
 837 //________________________________________________________________________
 838 void  AliAnalysisTaskEMCALPhoton::GetMcParts()
 839 {
 840    // Get MC particles.
 841   if(this->fDebug)
 842     printf("::GetMcParts() starting\n");
 843
 844   if (!this->fStack && !fAODMCParticles)
 845     return;
 846   if(this->fDebug)
 847     printf("either stack or aodmcpaticles exists\n");
 848   const AliVVertex *evtVtx = 0;
 849   if(this->fStack)
 850      evtVtx = fMCEvent->GetPrimaryVertex();
 851   else
 852     printf("no such thing as mc vvtx\n");
 853   /*if (!evtVtx)
 854     return;*/
 855   if(this->fDebug)
 856     printf("mc vvertex ok\n");
 857   Int_t nTracks = 0;
 858   if(this->fStack)
 859     nTracks = this->fStack->GetNtrack();
 860   else if(fAODMCParticles)
 861     nTracks = fAODMCParticles->GetEntriesFast();
 862   TParticle *mcP = 0;
 863   AliAODMCParticle *amcP = 0;
 864   if(this->fDebug)
 865     printf("loop in the %d mc particles starting\n",nTracks);
 866   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
 867     if(this->fStack)
 868       mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack));
 869     if(fAODMCParticles)
 870       amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
 871
 872     // primary particle
 873     Double_t dR = 0;
 874     if(mcP){
 875       dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
 876                        (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
 877                        (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
 878     }
 879     if((dR > 0.5))
 880       continue;
 881
 882
 883     // kinematic cuts
 884     Double_t pt = 0;
 885     Double_t eta = 0;
 886     Double_t phi  = 0;
 887     Int_t mother = -1;
 888     Int_t pdgcode = 0;
 889      if(mcP){
 890       pt = mcP->Pt() ;
 891       eta = mcP->Eta();
 892       phi = mcP->Phi();
 893       mother = mcP->GetMother(0);
 894       pdgcode = mcP->GetPdgCode();
 895      } else {
 896        if(amcP){
 897          pt = amcP->Pt();
 898          eta = amcP->Eta();
 899          phi = amcP->Phi();
 900          mother = amcP->GetMother();
 901          pdgcode = amcP->GetPdgCode();
 902        } else
 903        continue;
 904      }
 905     if (pt<0.5)
 906       continue;
 907
 908     if (TMath::Abs(eta)>0.7)
 909       continue;
 910
 911     if (phi<1.0||phi>3.3)
 912       continue;
 913
 914     if (mother!=6 && mother!=7 )
 915       continue;
 916
 917
 918     // pion or eta meson or direct photon
 919     if(pdgcode == 111) {
 920     } else if(pdgcode == 221) {
 921     } else if(pdgcode == 22 ) {
 922     } else {
 923       continue;
 924     }
 925
 926     FillMcPart( fMyMcIndex++, iTrack);
 927   }
 928   if(this->fDebug)
 929     printf("::GetMcParts() returning...\n\n");
 930 }
 931
 932 //________________________________________________________________________
 933 void  AliAnalysisTaskEMCALPhoton::FillMcPart(  Int_t itrack, Int_t label)
 934 {
 935   // Fill MC particles.
 936   Int_t nTracks = 0;
 937   TParticle *mcP = 0;
 938   AliAODMCParticle *amcP= 0;
 939   if(this->fStack){
 940     nTracks = this->fStack->GetNtrack();
 941     mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
 942   }
 943   else if(fAODMCParticles){
 944     nTracks = fAODMCParticles->GetEntriesFast();
 945     amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack));
 946  }
 947   if(this->fDebug)
 948     printf("\t::FillMcParts() starting with label %d\n", itrack);
 949    TVector3 vmcv;
 950    if(mcP)
 951      vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz());
 952    else {
 953      if(amcP)
 954        vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv());
 955      else
 956        return;
 957    }
 958
 959   AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
 960   mcp->fLabel = label ;
 961   if(mcP){
 962     mcp->fPdg = mcP->GetPdgCode() ;
 963     mcp->fPt = mcP->Pt() ;
 964     mcp->fEta = mcP->Eta() ;
 965     mcp->fPhi = mcP->Phi() ;
 966     mcp->fMother = mcP->GetMother(0) ;
 967     mcp->fFirstD = mcP->GetFirstDaughter() ;
 968     mcp->fLastD = mcP->GetLastDaughter() ;
 969     mcp->fStatus = mcP->GetStatusCode();
 970   }
 971   if(amcP){
 972     mcp->fPdg = amcP->GetPdgCode() ;
 973     mcp->fPt = amcP->Pt() ;
 974     mcp->fEta = amcP->Eta() ;
 975     mcp->fPhi = amcP->Phi() ;
 976     mcp->fMother = amcP->GetMother() ;
 977     mcp->fFirstD = amcP->GetDaughter(0) ;
 978     mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ;
 979     mcp->fStatus = amcP->GetStatus();
 980   }
 981   mcp->fVR = vmcv.Perp();
 982   if(vmcv.Perp()>0){
 983     mcp->fVEta = vmcv.Eta() ;
 984     mcp->fVPhi = vmcv.Phi() ;
 985   }
 986   if(itrack == 8){
 987     mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2);
 988     mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2);
 989     }
 990   if(this->fDebug)
 991     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);
 992   for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){
 993     if(id<=mcp->fMother)
 994       continue;
 995     if(id<0 || id>nTracks)
 996       continue;
 997     FillMcPart( fMyMcIndex++, id);
 998   }
 999
1000 }
1001 //________________________________________________________________________
1002 Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const
1003 {
1004   if(this->fDebug){
1005     printf("\t\t::GetMcIsolation() starting\n");
1006     //printf("\t\t   incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack);
1007   }
1008   if (!this->fStack && !this->fAODMCParticles){
1009     printf("\t\t\tNo MC stack/array!\n");
1010     return -1;
1011   }
1012   if(itrack<6 || itrack>8)
1013     return -1;
1014   if(this->fDebug)
1015     printf("\t\t\tparticle of interest selected\n");
1016   TParticle *mcP = 0;
1017   AliAODMCParticle *amcP = 0;
1018   Int_t pdgcode = 0;
1019   Float_t eta = 0;
1020   Float_t phi = 0;
1021   Double_t sumpt=0;
1022   Float_t dR;
1023   Int_t nparts =  0;
1024   if(this->fStack){
1025     nparts = fStack->GetNtrack();
1026     mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
1027     eta = mcP->Eta();
1028     phi = mcP->Phi();
1029     pdgcode = mcP->GetPdgCode();
1030   }
1031   if(this->fAODMCParticles){
1032     nparts = fAODMCParticles->GetEntriesFast();
1033     amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack));
1034     if(amcP){
1035       eta = amcP->Eta();
1036       phi = amcP->Phi();
1037       pdgcode = amcP->GetPdgCode();
1038     }
1039   }
1040   if(pdgcode!=22)
1041     return -1;
1042   TParticle *mcisop = 0;
1043   AliAODMCParticle *amcisop = 0;
1044   for(Int_t ip = 0; ip<nparts; ip++){
1045     if(ip==itrack)
1046       continue;
1047     if(this->fStack)
1048       mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip));
1049     if(fAODMCParticles)
1050       amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip));
1051     Int_t status = 0;
1052     Int_t mother = 0;
1053     Float_t pt = 0;
1054     Float_t isophi = -99;
1055     Float_t isoeta = -99;
1056     TVector3 vmcv;
1057     if(mcisop){
1058       status = mcisop->GetStatusCode();
1059       mother = mcisop->GetMother(0);
1060       pt = mcisop->Pt();
1061       isophi = mcisop->Phi();
1062       isoeta = mcisop->Eta();
1063       vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
1064     }
1065     else {
1066       if(amcisop){
1067         status = amcisop->GetStatus();
1068         mother = amcisop->GetMother();
1069         pt = amcisop->Pt();
1070         isophi = amcisop->Phi();
1071         isoeta = amcisop->Eta();
1072         vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv());
1073       }
1074       else
1075         continue;
1076     }
1077     if(status!=1)
1078       continue;
1079     if(mother == itrack)
1080       continue;
1081     if(pt<ptcut)
1082       continue;
1083     if(vmcv.Perp()>1)
1084       continue;
1085     dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta));
1086     if(dR>radius)
1087       continue;
1088     sumpt += pt;
1089   }
1090   if(this->fDebug)
1091     printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt);
1092   return sumpt;
1093  }
1094
1095 //________________________________________________________________________
1096 Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1097 {
1098   // Compute isolation based on tracks.
1099   if(this->fDebug)
1100     printf("\t::GetTrackIsolation() starting\n");
1101
1102   Double_t trkIsolation = 0;
1103   Double_t rad2 = radius*radius;
1104   Int_t ntrks = fSelPrimTracks->GetEntries();
1105   for(Int_t j = 0; j<ntrks; ++j) {
1106     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
1107     if (!track)
1108       continue;
1109     if (track->Pt()<pt)
1110       continue;
1111
1112     Float_t eta = track->Eta();
1113     Float_t phi = track->Phi();
1114     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1115     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1116     if(dist>rad2)
1117       continue;
1118     trkIsolation += track->Pt();
1119   }
1120   if(this->fDebug)
1121     printf("\t::GetTrackIsolation() returning\n\n");
1122   return trkIsolation;
1123 }
1124
1125 //________________________________________________________________________
1126 Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
1127 {
1128   // Get phi band.
1129
1130   if(!fSelPrimTracks)
1131     return 0;
1132   Int_t ntracks = fSelPrimTracks->GetEntries();
1133   Double_t loweta = eta - R;
1134   Double_t upeta = eta + R;
1135   Double_t upphi = phi + R;
1136   Double_t lowphi = phi - R;
1137   Double_t et = 0;
1138   for(Int_t itr=0; itr<ntracks; itr++){
1139     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr));
1140     if(!track)
1141       continue;
1142     if(track->Pt()<minpt)
1143       continue;
1144     if((track->Eta() < upeta) && (track->Eta() > loweta))
1145       continue;
1146     if((track->Phi() > upphi) || (track->Phi() < lowphi))
1147       continue;
1148     et+=track->Pt();
1149   }
1150   return et;
1151 }
1152
1153 //________________________________________________________________________
1154 Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1155 {
1156   // Calculate the energy of cross cells around the leading cell.
1157   if(!fVCells)
1158     return 0;
1159
1160   if (!fGeom)
1161     return 0;
1162
1163   Int_t iSupMod = -1;
1164   Int_t iTower  = -1;
1165   Int_t iIphi   = -1;
1166   Int_t iIeta   = -1;
1167   Int_t iphi    = -1;
1168   Int_t ieta    = -1;
1169   Int_t iphis   = -1;
1170   Int_t ietas   = -1;
1171
1172   Double_t crossEnergy = 0;
1173
1174   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1175   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1176
1177   Int_t ncells = cluster->GetNCells();
1178   for (Int_t i=0; i<ncells; i++) {
1179     Int_t cellAbsId = cluster->GetCellAbsId(i);
1180     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1181     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1182     Int_t aphidiff = TMath::Abs(iphi-iphis);
1183     if (aphidiff>1)
1184       continue;
1185     Int_t aetadiff = TMath::Abs(ieta-ietas);
1186     if (aetadiff>1)
1187       continue;
1188     if ( (aphidiff==1 && aetadiff==0) ||
1189         (aphidiff==0 && aetadiff==1) ) {
1190       crossEnergy += fVCells->GetCellAmplitude(cellAbsId);
1191     }
1192   }
1193
1194   return crossEnergy;
1195 }
1196
1197 //________________________________________________________________________
1198 Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1199 {
1200   // Get maximum energy of attached cell.
1201
1202   id = -1;
1203   if(!fVCells)
1204     return 0;
1205
1206   Double_t maxe = 0;
1207   Int_t ncells = cluster->GetNCells();
1208   for (Int_t i=0; i<ncells; i++) {
1209     Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1210     if (e>maxe) {
1211       maxe = e;
1212       id   = cluster->GetCellAbsId(i);
1213     }
1214   }
1215   return maxe;
1216 }
1217
1218 //________________________________________________________________________
1219 void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
1220 {
1221   // Called once at the end of the query
1222 /*  if(fIsGrid)
1223     return;*/
1224   if (fTree) {
1225     printf("***tree %s being saved***\n",fTree->GetName());
1226     TFile *f = OpenFile(2);
1227     TDirectory::TContext context(f);
1228     if (f)
1229       fTree->Write();
1230   }
1231 }