PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.cxx

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