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