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