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
  29
  30 #include "AliESDtrackCuts.h"
  31 #include "AliESDv0.h"
  32 #include "AliV0vertexer.h"
  33 #include "AliESDCaloCluster.h"
  34 #include "AliESDCaloCells.h"
  35 #include "AliEMCALGeometry.h"
  36 #include "AliEMCALRecoUtils.h"
  37 #include "AliVCluster.h"
  38 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  39 #include "TLorentzVector.h"
  40
  41
  42 #include "AliAnalysisTaskEMCALPhoton.h"
  43 #include "TFile.h"
  44
  45
  46 ClassImp(AliAnalysisTaskEMCALPhoton)
  47
  48 //________________________________________________________________________
  49 AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton(const char *name)
  50   : AliAnalysisTaskSE(name),
  51   fTrCuts(0),
  52   fPrTrCuts(0),
  53   fSelTracks(0),
  54   fSelPrimTracks(0),
  55   fPhotConvArray(0),
  56   fMyClusts(0),
  57   fMyAltClusts(0),
  58   fMyCells(0),
  59   fMyTracks(0),
  60   fMyMcParts(0),
  61   fHeader(0),
  62   fCaloClusters(0),
  63   fCaloClustersNew(0),
  64   fEMCalCells(0),
  65   fGeom(0x0),
  66   fTimeResTOF(0),
  67   fMipResponseTPC(0),
  68   fGeoName("EMCAL_COMPLETEV1"),
  69   fPeriod("LHC11c"),
  70   fIsTrain(0),
  71   fIsMC(0),
  72   fIsGrid(0),
  73   fClusThresh(2.),
  74   fClusterizer(0),
  75   fCaloClustersName("EmcalClusterv2"),
  76
  77
  78   fESD(0),
  79   fMCEvent(0),
  80   fStack(0x0),
  81
  82   fOutputList(0),
  83   fTree(0),
  84
  85
  86   fNV0sBefAndAftRerun(0),
  87   fConversionVtxXY(0),
  88   fInvMassV0(0),
  89   fInvMassV0KF(0),
  90   fInvMassV0SS(0),
  91   fDedxPAll(0)
  92
  93 {
  94   // Constructor
  95
  96   // Define input and output slots here
  97   // Input slot #0 works with a TChain
  98   DefineInput(0, TChain::Class());
  99   // Output slot #0 id reserved by the base class for AOD
 100   // Output slot #1 writes into a TH1 container
 101   DefineOutput(1, TList::Class());
 102 }
 103 //________________________________________________________________________
 104 void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects()
 105 {
 106   // Create histograms
 107   // Called once
 108
 109   fSelTracks = new TObjArray();
 110
 111   fSelPrimTracks = new TObjArray();
 112
 113   if (TClass::GetClass("AliPhotonHeaderObj"))
 114       TClass::GetClass("AliPhotonHeaderObj")->IgnoreTObjectStreamer();
 115   fHeader = new AliPhotonHeaderObj;
 116
 117   if (TClass::GetClass("AliPhotonConvObj"))
 118       TClass::GetClass("AliPhotonConvObj")->IgnoreTObjectStreamer();
 119   fPhotConvArray = new TClonesArray("AliPhotonConvObj");
 120
 121   if (TClass::GetClass("AliPhotonClusterObj"))
 122       TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
 123   fMyClusts = new TClonesArray("AliPhotonClusterObj");
 124
 125   if (TClass::GetClass("AliPhotonCellObj"))
 126       TClass::GetClass("AliPhotonCellObj")->IgnoreTObjectStreamer();
 127   fMyCells = new TClonesArray("AliPhotonCellObj");
 128
 129   if (TClass::GetClass("AliPhotonTrackObj"))
 130       TClass::GetClass("AliPhotonTrackObj")->IgnoreTObjectStreamer();
 131   fMyTracks = new TClonesArray("AliPhotonTrackObj");
 132
 133   if (fClusterizer || fIsTrain){
 134     if(fClusterizer)
 135       fCaloClustersName = fClusterizer->GetNewClusterArrayName();
 136     else {
 137       if(fPeriod.Contains("10h") || fPeriod.Contains("11h"))
 138         fCaloClustersName = "EmcalClustersv1";
 139       else
 140         fCaloClustersName = "EmcalClustersv2";
 141     }
 142     if (TClass::GetClass("AliPhotonClusterObj"))
 143         TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
 144     fMyAltClusts = new TClonesArray("AliPhotonClusterObj");
 145   }
 146   cout<<fCaloClustersName<<endl;
 147   if(fIsMC){
 148     if (TClass::GetClass("AliPhotonMcPartObj"))
 149         TClass::GetClass("AliPhotonMcPartObj")->IgnoreTObjectStreamer();
 150     fMyMcParts = new TClonesArray("AliPhotonMcPartObj");
 151   }
 152
 153   fCaloClusters = new TRefArray();
 154
 155   fOutputList = new TList();
 156   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 157
 158   if( !fTree){
 159     TFile *f = OpenFile(1);
 160     TDirectory::TContext context(f);
 161     if (f) {
 162       f->SetCompressionLevel(2);
 163       fTree = new TTree("photon_ana_out", "StandaloneTree");
 164       fTree->SetDirectory(f);
 165       if (fIsTrain) {
 166         fTree->SetAutoFlush(-2*1024*1024);
 167         fTree->SetAutoSave(0);
 168       } else {
 169         fTree->SetAutoFlush(-32*1024*1024);
 170         fTree->SetAutoSave(0);
 171       }
 172
 173     fTree->Branch("header", &fHeader, 16*1024, 99);
 174     fTree->Branch("conversions", &fPhotConvArray, 8*16*1024, 99);
 175     fTree->Branch("clusters", &fMyClusts, 8*16*1024, 99);
 176     if(fClusterizer || fIsTrain)
 177       fTree->Branch(fCaloClustersName, &fMyAltClusts, 8*16*1024, 99);
 178     fTree->Branch("cells", &fMyCells, 8*16*1024, 99);
 179     fTree->Branch("HighPtTracks", &fMyTracks, 8*16*1024, 99);
 180     if(fIsMC)
 181       fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99);
 182     }
 183   }
 184   if(fIsGrid)fOutputList->Add(fTree);
 185   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 186
 187
 188   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);
 189   fOutputList->Add(fNV0sBefAndAftRerun);
 190
 191   fConversionVtxXY = new TH2F("hConversionVtxXY","x and y of conversion vertex candidates;x;y",1000,-100,100,1000,-100,100);
 192   fOutputList->Add(fConversionVtxXY);
 193
 194   fInvMassV0 = new TH1F("hInvMassV0","v0->GetEffMass();v0->GetEffMass();dN/dM",400,0,4);
 195   fOutputList->Add(fInvMassV0);
 196
 197   fInvMassV0KF = new TH1F("hInvMassV0KF","Inv. mass calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
 198   fOutputList->Add(fInvMassV0KF);
 199
 200   fInvMassV0SS = new TH1F("hInvMassV0SS","Inv. mass (same sign) calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
 201   fOutputList->Add(fInvMassV0SS);
 202
 203   fDedxPAll = new TH2F("hDedxPAll","dE/dx vs p (all selected tracks);p (GeV/c);dE/dx (a.u.)",400,0,40, 150,0,150);
 204   fOutputList->Add(fDedxPAll);
 205
 206
 207   PostData(1, fOutputList);
 208 }
 209
 210 //________________________________________________________________________
 211 void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *)
 212 {
 213   //event trigger selection
 214   Bool_t isSelected = 0;
 215   if(fPeriod.Contains("11")){
 216     if(fPeriod.Contains("11a"))
 217       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 218     if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") )
 219       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 220     if(fPeriod.Contains("11h") )
 221       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA);
 222
 223   }
 224   else
 225     isSelected = 1;
 226 /*  if(!isSelected )
 227         return;   */
 228   // Main loop
 229   // Called for each event
 230
 231   // Post output data.
 232   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 233   if (!fESD) {
 234     printf("ERROR: fESD not available, returning...\n");
 235     return;
 236   }
 237
 238
 239
 240   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 241   if(!pv)
 242     return;
 243   if(TMath::Abs(pv->GetZ())>15)
 244     return;
 245   //fCaloClustersNew = dynamic_cast<TClonesArray*>(fESD->FindListObject(fCaloClustersName));
 246
 247
 248   // Track loop to fill a pT spectrum
 249   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 250     AliESDtrack* track = fESD->GetTrack(iTracks);
 251     if (!track)
 252       continue;
 253
 254
 255     if (fTrCuts && fTrCuts->IsSelected(track)){
 256       fSelTracks->Add(track);
 257       fDedxPAll->Fill(track->P(), track->GetTPCsignal());
 258     }
 259     if (fPrTrCuts && fPrTrCuts->IsSelected(track))
 260       fSelPrimTracks->Add(track);
 261
 262
 263
 264   } //track loop
 265   fHeader->fTrClassMask    = fESD->GetHeader()->GetTriggerMask();
 266   fHeader->fTrCluster      = fESD->GetHeader()->GetTriggerCluster();
 267   AliCentrality *cent = InputEvent()->GetCentrality();
 268   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
 269   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
 270   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
 271   if(!fIsTrain){
 272     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 273       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 274         break;
 275       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 276     }
 277   }
 278   fESD->GetEMCALClusters(fCaloClusters);
 279   fHeader->fNClus = fCaloClusters->GetEntries();
 280   fEMCalCells = fESD->GetEMCALCells();
 281   fHeader->fNCells = fEMCalCells->GetNumberOfCells();
 282   fMCEvent = MCEvent();
 283   if(fMCEvent)
 284     fStack = (AliStack*)fMCEvent->Stack();
 285
 286
 287   FindConversions();
 288 //   cout<<"****calling FillMyCells()...\n";
 289   FillMyCells();
 290 //   cout<<"****calling FillMyClus...\n";
 291   FillMyClusters();
 292   if(fCaloClustersNew)
 293     FillMyAltClusters();
 294 //   cout<<"****calling FillHighPt...\n";
 295   FillHighPtTracks();
 296   if(fIsMC)
 297     GetMcParts();
 298
 299   fTree->Fill();
 300   fSelTracks->Clear();
 301   fSelPrimTracks->Clear();
 302   fPhotConvArray->Clear();
 303   fMyClusts->Clear();
 304   if(fMyAltClusts)
 305     fMyAltClusts->Clear();
 306   fMyCells->Clear();
 307   fMyTracks->Clear();
 308   if(fMyMcParts)
 309     fMyMcParts->Clear();
 310   fCaloClusters->Clear();
 311   if(fCaloClustersNew)
 312     fCaloClustersNew->Clear();
 313   PostData(1, fOutputList);
 314 }
 315
 316 //________________________________________________________________________
 317 void AliAnalysisTaskEMCALPhoton::FindConversions()
 318 {
 319   if(!fTrCuts)
 320     return;
 321   Int_t iconv = 0;
 322   Int_t nV0Orig = fESD->GetNumberOfV0s();
 323   Int_t ntracks = fESD->GetNumberOfTracks();
 324   fESD->ResetV0s();
 325   AliV0vertexer lV0vtxer;
 326   lV0vtxer.Tracks2V0vertices(fESD);
 327   Int_t nV0s = fESD->GetNumberOfV0s();
 328   fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
 329   for(Int_t iv0=0; iv0<nV0s; iv0++){
 330     AliESDv0 * v0 = fESD->GetV0(iv0);
 331     if(!v0)
 332       continue;
 333     Double_t vpos[3];
 334     fInvMassV0->Fill(v0->GetEffMass());
 335     v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
 336     Int_t ipos = v0->GetPindex();
 337     Int_t ineg = v0->GetNindex();
 338     if(ipos<0 || ipos> ntracks)
 339       continue;
 340     if(ineg<0 || ineg> ntracks)
 341       continue;
 342     AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
 343     AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
 344     if(!pos || !neg)
 345       continue;
 346     if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
 347       continue;
 348     /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
 349       continue;*/
 350     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
 351     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
 352     if(!paramPos || !paramNeg)
 353       continue;
 354     if(pos->GetSign() <0){//change tracks
 355       pos=neg ;
 356       neg=fESD->GetTrack(v0->GetPindex()) ;
 357       paramPos=paramNeg ;
 358       paramNeg=v0->GetParamP() ;
 359     }
 360     AliKFParticle negKF(*paramNeg,11);
 361     AliKFParticle posKF(*paramPos,-11);
 362     AliKFParticle photKF(negKF,posKF) ;
 363
 364
 365     if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
 366       continue ;
 367
 368
 369     if(pos->GetSign() == neg->GetSign()){
 370       fInvMassV0SS->Fill(photKF.GetMass());
 371       continue;
 372     }
 373     fInvMassV0KF->Fill(photKF.GetMass());
 374     TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
 375     if(photLV.M()>0.05 || photLV.M()<0)
 376       continue;
 377     fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
 378     Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
 379     if(convPhi<0)
 380       convPhi+=TMath::TwoPi();
 381     TVector3 vecpos(vpos);
 382     Double_t v0Phi = 0;
 383     if(vecpos.Perp()>0)
 384       vecpos.Phi();
 385     if(v0Phi<0)
 386       v0Phi+=TMath::TwoPi();
 387     AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
 388     myconv->fPt = photLV.Pt();
 389     myconv->fEta = photLV.Eta();
 390     myconv->fPhi = convPhi;
 391     myconv->fVR = vecpos.Perp();
 392     if(vecpos.Perp()>0)
 393       myconv->fVEta = vecpos.Eta();
 394     myconv->fVPhi = v0Phi;
 395     myconv->fMass = photLV.M();
 396     myconv->fMcLabel = -3; //WARNING: include the correct labeling
 397     //negative daughter
 398    myconv->fNegPt      =  negKF.GetPt();
 399    myconv->fNegEta     =  negKF.GetEta();
 400    Double_t trackPhi   =  negKF.GetPhi();
 401    if(trackPhi<0)
 402      trackPhi+=TMath::TwoPi();
 403    myconv->fNegPhi     =  trackPhi;
 404    myconv->fNegDedx    =  neg->GetTPCsignal();
 405    myconv->fNegMcLabel =  neg->GetLabel();
 406     //negative daughter
 407    myconv->fPosPt      =  posKF.GetPt();
 408    myconv->fPosEta     =  posKF.GetEta();
 409    trackPhi            =  posKF.GetPhi();
 410    if(trackPhi<0)
 411      trackPhi+=TMath::TwoPi();
 412    myconv->fPosPhi     =  trackPhi;
 413    myconv->fPosDedx    =  pos->GetTPCsignal();
 414    myconv->fPosMcLabel =  pos->GetLabel();
 415   }
 416   return;
 417 }
 418 //________________________________________________________________________
 419 void AliAnalysisTaskEMCALPhoton::FillMyCells()
 420 {
 421   if(!fEMCalCells)
 422     return;
 423 //   cout<<"++++FillMyCells starts!\n";
 424   Int_t ncells = fEMCalCells->GetNumberOfCells();
 425   Int_t mcel = 0;
 426   for(Int_t icell = 0; icell<ncells; icell++){
 427     Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 428     AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
 429     Float_t eta=-1, phi=-1;
 430     if(!fGeom){
 431       cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
 432       return;
 433     }
 434     if(!fGeom)
 435       return;
 436     if(!fIsMC)fGeom->EtaPhiFromIndex(absID,eta,phi);
 437     Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 438     mycell->fAbsID = absID;
 439     mycell->fE = fEMCalCells->GetCellAmplitude(absID);
 440     mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
 441     mycell->fEta = eta;
 442     mycell->fPhi = phi;
 443     mycell->fTime = fEMCalCells->GetCellTime(absID);
 444   }
 445 //   cout<<"++++FillMyCells done!\n";
 446 }
 447 //________________________________________________________________________
 448 void AliAnalysisTaskEMCALPhoton::FillMyClusters()
 449 {
 450   if(!fCaloClusters)
 451     return;
 452   Int_t nclus = fCaloClusters->GetEntries();
 453   Int_t mcl = 0;
 454   for(Int_t ic=0; ic < nclus; ic++){
 455     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 456     if(!clus)
 457       continue;
 458     if(!clus->IsEMCAL())
 459       continue;
 460     if(clus->E() < fClusThresh)
 461       continue;
 462     Float_t pos[3];
 463     clus->GetPosition(pos);
 464     TVector3 cpos(pos);
 465     TString cellsAbsID;
 466     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
 467     myclus->fE       = clus->E();
 468     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 469     myclus->fR       = cpos.Perp();
 470     myclus->fEta     = cpos.Eta();
 471     myclus->fPhi     = cpos.Phi();
 472     if(cpos.Phi()<0){
 473       myclus->fPhi+=TMath::TwoPi();
 474     }
 475     myclus->fN       = clus->GetNCells();
 476     Short_t  id = -1;
 477     myclus->fEmax    = GetMaxCellEnergy( clus, id);
 478     myclus->fIdmax   = id;
 479     myclus->fTmax    = fEMCalCells->GetCellTime(id);
 480     myclus->fEcross  = GetCrossEnergy( clus, id);
 481     myclus->fDisp    = clus->GetDispersion();
 482     myclus->fM20     = clus->GetM20();
 483     myclus->fM02     = clus->GetM02();
 484     myclus->fTrDEta  = clus->GetTrackDz();
 485     myclus->fTrDPhi  = clus->GetTrackDx();
 486     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 487     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 488     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 489     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 490     myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 491     myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 492     myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 493     myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 494     for(Int_t icell=0;icell<myclus->fN;icell++){
 495       Int_t absID = clus->GetCellAbsId(icell);
 496       cellsAbsID.Append(Form("%d",absID));
 497       cellsAbsID.Append(";");
 498     }
 499     myclus->fCellsAbsId = cellsAbsID;
 500     myclus->fMcLabel = clus->GetLabel();
 501     Int_t matchIndex = clus->GetTrackMatchedIndex();
 502     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
 503       myclus->fTrEp = -1;
 504       continue;
 505     }
 506     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
 507     if(!track){
 508       myclus->fTrEp = -1;
 509       continue;
 510     }
 511     if(!fPrTrCuts){
 512       myclus->fTrEp = -1;
 513       continue;
 514     }
 515     if(!fPrTrCuts->IsSelected(track)){
 516       myclus->fTrEp = -1;
 517       continue;
 518     }
 519     myclus->fTrEp = clus->E()/track->P();
 520   }
 521
 522 }
 523 //________________________________________________________________________
 524 void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
 525 {
 526   if(!fCaloClustersNew)
 527     return;
 528   Int_t nclus = fCaloClustersNew->GetEntries();
 529   Int_t mcl = 0;
 530   for(Int_t ic=0; ic < nclus; ic++){
 531     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
 532     if(!clus)
 533       continue;
 534     if(!clus->IsEMCAL())
 535       continue;
 536     if(clus->E() < fClusThresh)
 537       continue;
 538     Float_t pos[3];
 539     clus->GetPosition(pos);
 540     TVector3 cpos(pos);
 541     TString cellsAbsID;
 542     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
 543     myclus->fE       = clus->E();
 544     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 545     myclus->fR       = cpos.Perp();
 546     myclus->fEta     = cpos.Eta();
 547     myclus->fPhi     = cpos.Phi();
 548     if(cpos.Phi()<0){
 549       myclus->fPhi+=TMath::TwoPi();
 550     }
 551     myclus->fN       = clus->GetNCells();
 552     myclus->fDisp    = clus->GetDispersion();
 553     myclus->fM20     = clus->GetM20();
 554     myclus->fM02     = clus->GetM02();
 555     myclus->fTrDEta  = clus->GetTrackDz();
 556     myclus->fTrDPhi  = clus->GetTrackDx();
 557     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 558     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 559     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 560     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 561     for(Int_t icell=0;icell<myclus->fN;icell++){
 562       Int_t absID = clus->GetCellAbsId(icell);
 563       cellsAbsID.Append(Form("%d",absID));
 564       cellsAbsID.Append(";");
 565     }
 566     myclus->fCellsAbsId = cellsAbsID;
 567     myclus->fMcLabel = clus->GetLabel();
 568     Int_t matchIndex = clus->GetTrackMatchedIndex();
 569     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
 570       myclus->fTrEp = -1;
 571       continue;
 572     }
 573     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
 574     if(!track){
 575       myclus->fTrEp = -1;
 576       continue;
 577     }
 578     if(!fPrTrCuts){
 579       myclus->fTrEp = -1;
 580       continue;
 581     }
 582     if(!fPrTrCuts->IsSelected(track)){
 583       myclus->fTrEp = -1;
 584       continue;
 585     }
 586     myclus->fTrEp = clus->E()/track->P();
 587   }
 588
 589 }
 590 //________________________________________________________________________
 591 void  AliAnalysisTaskEMCALPhoton::FillHighPtTracks()
 592 {
 593   if(!fSelPrimTracks)
 594     return;
 595   Int_t ntracks = fSelPrimTracks->GetEntries();
 596   Int_t imtr = 0;
 597   for(Int_t it=0;it<ntracks; it++){
 598     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
 599     if(!track)
 600       continue;
 601     if(track->Pt()<3)
 602       continue;
 603     AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
 604     mtr->fPt = track->Pt();
 605     mtr->fEta = track->Eta();
 606     mtr->fPhi = track->Phi();
 607     mtr->fDedx = track->GetTPCsignal();
 608     mtr->fMcLabel = track->GetLabel();
 609   }
 610 }
 611 //________________________________________________________________________
 612 void  AliAnalysisTaskEMCALPhoton::GetMcParts()
 613 {
 614   if (!fStack)
 615     return;
 616 //   cout<<"++++Get - MC starts!\n";
 617
 618   const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
 619   if (!evtVtx)
 620     return;
 621   Int_t ipart = 0;
 622   Int_t nTracks = fStack->GetNtrack();
 623   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
 624     TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
 625     if (!mcP)
 626       continue;
 627     // primary particle
 628     Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
 629                               (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
 630                               (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
 631     if(dR > 0.5)
 632       continue;
 633
 634     // kinematic cuts
 635     Double_t pt = mcP->Pt() ;
 636     if (pt<0.5)
 637       continue;
 638     Double_t eta = mcP->Eta();
 639     if (eta<-0.7||eta>0.7)
 640       continue;
 641     Double_t phi  = mcP->Phi();
 642     if (phi<1.0||phi>3.3)
 643       continue;
 644     // pion or eta meson or direct photon
 645     if(mcP->GetPdgCode() == 111) {
 646     } else if(mcP->GetPdgCode() == 221) {
 647     } else if(mcP->GetPdgCode() == 22 ) {
 648     } else
 649       continue;
 650     bool checkIfAlreadySaved = false;
 651     for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
 652       AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
 653       if(!mymc)
 654         continue;
 655       if(mymc->fLabel == iTrack)
 656         checkIfAlreadySaved = true;
 657     }
 658     if(!checkIfAlreadySaved)
 659       FillMcPart(mcP, ipart++, iTrack);
 660     for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
 661       if(id<=mcP->GetMother(0))
 662         continue;
 663       if(id<0 || id>nTracks)
 664         continue;
 665       TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
 666       if(!mcD)
 667         continue;
 668       FillMcPart(mcD, ipart++, id);
 669     }
 670   }
 671 //   cout<<"++++Get - MC done!\n";
 672
 673 }
 674 //________________________________________________________________________
 675 void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
 676 {
 677   if(!mcP)
 678     return;
 679   TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
 680   AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
 681   mcp->fLabel = itrack ;
 682   mcp->fPdg = mcP->GetPdgCode() ;
 683   mcp->fPt = mcP->Pt() ;
 684   mcp->fEta = mcP->Eta() ;
 685   mcp->fPhi = mcP->Phi() ;
 686   mcp->fVR = vmcv.Perp();
 687   if(vmcv.Perp()>0){
 688     mcp->fVEta = vmcv.Eta() ;
 689     mcp->fVPhi = vmcv.Phi() ;
 690   }
 691   mcp->fMother = mcP->GetMother(0) ;
 692 }
 693 //________________________________________________________________________
 694 Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
 695 {
 696   // Compute isolation based on tracks.
 697
 698   Double_t trkIsolation = 0;
 699   Double_t rad2 = radius*radius;
 700   Int_t ntrks = fSelPrimTracks->GetEntries();
 701   for(Int_t j = 0; j<ntrks; ++j) {
 702     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
 703     if (!track)
 704       continue;
 705     if (track->Pt()<pt)
 706       continue;
 707
 708     Float_t eta = track->Eta();
 709     Float_t phi = track->Phi();
 710     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
 711     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
 712     if(dist>rad2)
 713       continue;
 714     trkIsolation += track->Pt();
 715   }
 716   return trkIsolation;
 717 }
 718 //________________________________________________________________________
 719 Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
 720 {
 721   if(!fSelPrimTracks)
 722     return 0;
 723   Int_t ntracks = fSelPrimTracks->GetEntries();
 724   Double_t loweta = eta - R;
 725   Double_t upeta = eta + R;
 726   Double_t upphi = phi + R;
 727   Double_t lowphi = phi - R;
 728   Double_t et = 0;
 729   for(Int_t itr=0; itr<ntracks; itr++){
 730     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
 731     if(!track)
 732       continue;
 733     if(track->Pt()<minpt)
 734       continue;
 735     if((track->Eta() < upeta) && (track->Eta() > loweta))
 736       continue;
 737     if((track->Phi() > upphi) || (track->Phi() < lowphi))
 738       continue;
 739     et+=track->Pt();
 740   }
 741   return et;
 742 }
 743 //________________________________________________________________________
 744 Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 745 {
 746   // Calculate the energy of cross cells around the leading cell.
 747
 748   AliVCaloCells *cells = 0;
 749   cells = fESD->GetEMCALCells();
 750   if (!cells)
 751     return 0;
 752
 753
 754   if (!fGeom)
 755     return 0;
 756
 757   Int_t iSupMod = -1;
 758   Int_t iTower  = -1;
 759   Int_t iIphi   = -1;
 760   Int_t iIeta   = -1;
 761   Int_t iphi    = -1;
 762   Int_t ieta    = -1;
 763   Int_t iphis   = -1;
 764   Int_t ietas   = -1;
 765
 766   Double_t crossEnergy = 0;
 767
 768   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 769   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 770
 771   Int_t ncells = cluster->GetNCells();
 772   for (Int_t i=0; i<ncells; i++) {
 773     Int_t cellAbsId = cluster->GetCellAbsId(i);
 774     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 775     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 776     Int_t aphidiff = TMath::Abs(iphi-iphis);
 777     if (aphidiff>1)
 778       continue;
 779     Int_t aetadiff = TMath::Abs(ieta-ietas);
 780     if (aetadiff>1)
 781       continue;
 782     if ( (aphidiff==1 && aetadiff==0) ||
 783         (aphidiff==0 && aetadiff==1) ) {
 784       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 785     }
 786   }
 787
 788   return crossEnergy;
 789 }
 790
 791
 792
 793 //________________________________________________________________________
 794 Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 795 {
 796   // Get maximum energy of attached cell.
 797
 798   id = -1;
 799
 800   AliVCaloCells *cells = 0;
 801   cells = fESD->GetEMCALCells();
 802   if (!cells)
 803     return 0;
 804
 805   Double_t maxe = 0;
 806   Int_t ncells = cluster->GetNCells();
 807   for (Int_t i=0; i<ncells; i++) {
 808     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 809     if (e>maxe) {
 810       maxe = e;
 811       id   = cluster->GetCellAbsId(i);
 812     }
 813   }
 814   return maxe;
 815 }
 816 //________________________________________________________________________
 817 void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
 818 {
 819   // Draw result to the screen
 820   // Called once at the end of the query
 821     if(fIsGrid)
 822       return;
 823     if (fTree) {
 824     TFile *f = OpenFile(1);
 825     TDirectory::TContext context(f);
 826     if (f)
 827       fTree->Write();
 828   }
 829
 830
 831 }