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("HighPtTracks", &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   else
 266     isSelected = 1;
 267
 268   // Post output data.
 269   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 270   if (!fESD) {
 271     printf("ERROR: fESD not available, returning...\n");
 272     return;
 273   }
 274
 275   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 276   if(!pv)
 277     return;
 278   if(TMath::Abs(pv->GetZ())>15)
 279     return;
 280
 281   // Track loop to fill a pT spectrum
 282   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 283     AliESDtrack* track = fESD->GetTrack(iTracks);
 284     if (!track)
 285       continue;
 286
 287
 288     if (fTrCuts && fTrCuts->IsSelected(track)){
 289       fSelTracks->Add(track);
 290       fDedxPAll->Fill(track->P(), track->GetTPCsignal());
 291     }
 292     if (fPrTrCuts && fPrTrCuts->IsSelected(track))
 293       fSelPrimTracks->Add(track);
 294   } //track loop
 295
 296   fHeader->fTrClassMask    = fESD->GetHeader()->GetTriggerMask();
 297   fHeader->fTrCluster      = fESD->GetHeader()->GetTriggerCluster();
 298   AliCentrality *cent = InputEvent()->GetCentrality();
 299   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
 300   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
 301   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
 302   if(!fIsTrain){
 303     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 304       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 305         break;
 306       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 307     }
 308   }
 309   fESD->GetEMCALClusters(fCaloClusters);
 310   fHeader->fNClus = fCaloClusters->GetEntries();
 311   fEMCalCells = fESD->GetEMCALCells();
 312   fHeader->fNCells = fEMCalCells->GetNumberOfCells();
 313   fMCEvent = MCEvent();
 314   if(fMCEvent)
 315     fStack = (AliStack*)fMCEvent->Stack();
 316
 317
 318   FindConversions();
 319   FillMyCells();
 320   FillMyClusters();
 321   if(fCaloClustersNew)
 322     FillMyAltClusters();
 323   FillHighPtTracks();
 324   if(fIsMC)
 325     GetMcParts();
 326
 327   fTree->Fill();
 328   fSelTracks->Clear();
 329   fSelPrimTracks->Clear();
 330   fPhotConvArray->Clear();
 331   fMyClusts->Clear();
 332   if(fMyAltClusts)
 333     fMyAltClusts->Clear();
 334   fMyCells->Clear();
 335   fMyTracks->Clear();
 336   if(fMyMcParts)
 337     fMyMcParts->Clear();
 338   fCaloClusters->Clear();
 339   if(fCaloClustersNew)
 340     fCaloClustersNew->Clear();
 341   PostData(1, fOutputList);
 342 }
 343
 344 //________________________________________________________________________
 345 void AliAnalysisTaskEMCALPhoton::FindConversions()
 346 {
 347   // Find conversion.
 348
 349   if(!fTrCuts)
 350     return;
 351   Int_t iconv = 0;
 352   Int_t nV0Orig = fESD->GetNumberOfV0s();
 353   Int_t ntracks = fESD->GetNumberOfTracks();
 354   fESD->ResetV0s();
 355   AliV0vertexer lV0vtxer;
 356   lV0vtxer.Tracks2V0vertices(fESD);
 357   Int_t nV0s = fESD->GetNumberOfV0s();
 358   fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
 359   for(Int_t iv0=0; iv0<nV0s; iv0++){
 360     AliESDv0 * v0 = fESD->GetV0(iv0);
 361     if(!v0)
 362       continue;
 363     Double_t vpos[3];
 364     fInvMassV0->Fill(v0->GetEffMass());
 365     v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
 366     Int_t ipos = v0->GetPindex();
 367     Int_t ineg = v0->GetNindex();
 368     if(ipos<0 || ipos> ntracks)
 369       continue;
 370     if(ineg<0 || ineg> ntracks)
 371       continue;
 372     AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
 373     AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
 374     if(!pos || !neg)
 375       continue;
 376     if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
 377       continue;
 378     /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
 379       continue;*/
 380     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
 381     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
 382     if(!paramPos || !paramNeg)
 383       continue;
 384     if(pos->GetSign() <0){//change tracks
 385       pos=neg ;
 386       neg=fESD->GetTrack(v0->GetPindex()) ;
 387       paramPos=paramNeg ;
 388       paramNeg=v0->GetParamP() ;
 389     }
 390     AliKFParticle negKF(*paramNeg,11);
 391     AliKFParticle posKF(*paramPos,-11);
 392     AliKFParticle photKF(negKF,posKF) ;
 393
 394     if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
 395       continue ;
 396
 397     if(pos->GetSign() == neg->GetSign()){
 398       fInvMassV0SS->Fill(photKF.GetMass());
 399       continue;
 400     }
 401     fInvMassV0KF->Fill(photKF.GetMass());
 402     TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
 403     if(photLV.M()>0.05 || photLV.M()<0)
 404       continue;
 405     fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
 406     Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
 407     if(convPhi<0)
 408       convPhi+=TMath::TwoPi();
 409     TVector3 vecpos(vpos);
 410     Double_t v0Phi = 0;
 411     if(vecpos.Perp()>0)
 412       vecpos.Phi();
 413     if(v0Phi<0)
 414       v0Phi+=TMath::TwoPi();
 415     AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
 416     myconv->fPt = photLV.Pt();
 417     myconv->fEta = photLV.Eta();
 418     myconv->fPhi = convPhi;
 419     myconv->fVR = vecpos.Perp();
 420     if(vecpos.Perp()>0)
 421       myconv->fVEta = vecpos.Eta();
 422     myconv->fVPhi = v0Phi;
 423     myconv->fMass = photLV.M();
 424     myconv->fMcLabel = -3; //WARNING: include the correct labeling
 425     //negative daughter
 426    myconv->fNegPt      =  negKF.GetPt();
 427    myconv->fNegEta     =  negKF.GetEta();
 428    Double_t trackPhi   =  negKF.GetPhi();
 429    if(trackPhi<0)
 430      trackPhi+=TMath::TwoPi();
 431    myconv->fNegPhi     =  trackPhi;
 432    myconv->fNegDedx    =  neg->GetTPCsignal();
 433    myconv->fNegMcLabel =  neg->GetLabel();
 434     //negative daughter
 435    myconv->fPosPt      =  posKF.GetPt();
 436    myconv->fPosEta     =  posKF.GetEta();
 437    trackPhi            =  posKF.GetPhi();
 438    if(trackPhi<0)
 439      trackPhi+=TMath::TwoPi();
 440    myconv->fPosPhi     =  trackPhi;
 441    myconv->fPosDedx    =  pos->GetTPCsignal();
 442    myconv->fPosMcLabel =  pos->GetLabel();
 443   }
 444   return;
 445 }
 446
 447 //________________________________________________________________________
 448 void AliAnalysisTaskEMCALPhoton::FillMyCells()
 449 {
 450   // Fill cells.
 451
 452   if (!fEMCalCells)
 453     return;
 454   Int_t ncells = fEMCalCells->GetNumberOfCells();
 455   Int_t mcel = 0;
 456   for(Int_t icell = 0; icell<ncells; icell++){
 457     Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 458     AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
 459     Float_t eta=-1, phi=-1;
 460     if(!fGeom){
 461       cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
 462       return;
 463     }
 464     if(!fGeom)
 465       return;
 466     if(!fIsMC)fGeom->EtaPhiFromIndex(absID,eta,phi);
 467     Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 468     mycell->fAbsID = absID;
 469     mycell->fE = fEMCalCells->GetCellAmplitude(absID);
 470     mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
 471     mycell->fEta = eta;
 472     mycell->fPhi = phi;
 473     mycell->fTime = fEMCalCells->GetCellTime(absID);
 474   }
 475 }
 476
 477 //________________________________________________________________________
 478 void AliAnalysisTaskEMCALPhoton::FillMyClusters()
 479 {
 480   // Fill clusters.
 481
 482   if (!fCaloClusters)
 483     return;
 484   Int_t nclus = fCaloClusters->GetEntries();
 485   Int_t mcl = 0;
 486   for(Int_t ic=0; ic < nclus; ic++){
 487     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 488     if(!clus)
 489       continue;
 490     if(!clus->IsEMCAL())
 491       continue;
 492     if(clus->E() < fClusThresh)
 493       continue;
 494     Float_t pos[3];
 495     clus->GetPosition(pos);
 496     TVector3 cpos(pos);
 497     TString cellsAbsID;
 498     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
 499     myclus->fE       = clus->E();
 500     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 501     myclus->fR       = cpos.Perp();
 502     myclus->fEta     = cpos.Eta();
 503     myclus->fPhi     = cpos.Phi();
 504     if(cpos.Phi()<0){
 505       myclus->fPhi+=TMath::TwoPi();
 506     }
 507     myclus->fN       = clus->GetNCells();
 508     Short_t  id = -1;
 509     myclus->fEmax    = GetMaxCellEnergy( clus, id);
 510     myclus->fIdmax   = id;
 511     myclus->fTmax    = fEMCalCells->GetCellTime(id);
 512     myclus->fEcross  = GetCrossEnergy( clus, id);
 513     myclus->fDisp    = clus->GetDispersion();
 514     myclus->fM20     = clus->GetM20();
 515     myclus->fM02     = clus->GetM02();
 516     myclus->fTrDEta  = clus->GetTrackDz();
 517     myclus->fTrDPhi  = clus->GetTrackDx();
 518     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 519     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 520     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 521     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 522     myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 523     myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 524     myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 525     myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 526     for(Int_t icell=0;icell<myclus->fN;icell++){
 527       Int_t absID = clus->GetCellAbsId(icell);
 528       cellsAbsID.Append(Form("%d",absID));
 529       cellsAbsID.Append(";");
 530     }
 531     myclus->fCellsAbsId = cellsAbsID;
 532     myclus->fMcLabel = clus->GetLabel();
 533     Int_t matchIndex = clus->GetTrackMatchedIndex();
 534     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
 535       myclus->fTrEp = -1;
 536       continue;
 537     }
 538     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
 539     if(!track){
 540       myclus->fTrEp = -1;
 541       continue;
 542     }
 543     if(!fPrTrCuts){
 544       myclus->fTrEp = -1;
 545       continue;
 546     }
 547     if(!fPrTrCuts->IsSelected(track)){
 548       myclus->fTrEp = -1;
 549       continue;
 550     }
 551     myclus->fTrEp = clus->E()/track->P();
 552   }
 553
 554 }
 555 //________________________________________________________________________
 556 void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
 557 {
 558   // Fill clusters.
 559
 560   if(!fCaloClustersNew)
 561     return;
 562   Int_t nclus = fCaloClustersNew->GetEntries();
 563   Int_t mcl = 0;
 564   for(Int_t ic=0; ic < nclus; ic++){
 565     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
 566     if(!clus)
 567       continue;
 568     if(!clus->IsEMCAL())
 569       continue;
 570     if(clus->E() < fClusThresh)
 571       continue;
 572     Float_t pos[3];
 573     clus->GetPosition(pos);
 574     TVector3 cpos(pos);
 575     TString cellsAbsID;
 576     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
 577     myclus->fE       = clus->E();
 578     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
 579     myclus->fR       = cpos.Perp();
 580     myclus->fEta     = cpos.Eta();
 581     myclus->fPhi     = cpos.Phi();
 582     if(cpos.Phi()<0){
 583       myclus->fPhi+=TMath::TwoPi();
 584     }
 585     myclus->fN       = clus->GetNCells();
 586     myclus->fDisp    = clus->GetDispersion();
 587     myclus->fM20     = clus->GetM20();
 588     myclus->fM02     = clus->GetM02();
 589     myclus->fTrDEta  = clus->GetTrackDz();
 590     myclus->fTrDPhi  = clus->GetTrackDx();
 591     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
 592     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
 593     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
 594     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
 595     for(Int_t icell=0;icell<myclus->fN;icell++){
 596       Int_t absID = clus->GetCellAbsId(icell);
 597       cellsAbsID.Append(Form("%d",absID));
 598       cellsAbsID.Append(";");
 599     }
 600     myclus->fCellsAbsId = cellsAbsID;
 601     myclus->fMcLabel = clus->GetLabel();
 602     Int_t matchIndex = clus->GetTrackMatchedIndex();
 603     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
 604       myclus->fTrEp = -1;
 605       continue;
 606     }
 607     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
 608     if(!track){
 609       myclus->fTrEp = -1;
 610       continue;
 611     }
 612     if(!fPrTrCuts){
 613       myclus->fTrEp = -1;
 614       continue;
 615     }
 616     if(!fPrTrCuts->IsSelected(track)){
 617       myclus->fTrEp = -1;
 618       continue;
 619     }
 620     myclus->fTrEp = clus->E()/track->P();
 621   }
 622
 623 }
 624 //________________________________________________________________________
 625 void  AliAnalysisTaskEMCALPhoton::FillHighPtTracks()
 626 {
 627   // Fill high pt tracks.
 628
 629   if(!fSelPrimTracks)
 630     return;
 631   Int_t ntracks = fSelPrimTracks->GetEntries();
 632   Int_t imtr = 0;
 633   for(Int_t it=0;it<ntracks; it++){
 634     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
 635     if(!track)
 636       continue;
 637     if(track->Pt()<3)
 638       continue;
 639     AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
 640     mtr->fPt = track->Pt();
 641     mtr->fEta = track->Eta();
 642     mtr->fPhi = track->Phi();
 643     mtr->fDedx = track->GetTPCsignal();
 644     mtr->fMcLabel = track->GetLabel();
 645   }
 646 }
 647
 648 //________________________________________________________________________
 649 void  AliAnalysisTaskEMCALPhoton::GetMcParts()
 650 {
 651   // Get MC particles.
 652
 653   if (!fStack)
 654     return;
 655
 656   const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
 657   if (!evtVtx)
 658     return;
 659   Int_t ipart = 0;
 660   Int_t nTracks = fStack->GetNtrack();
 661   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
 662     TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
 663     if (!mcP)
 664       continue;
 665     // primary particle
 666     Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
 667                               (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
 668                               (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
 669     if(dR > 0.5)
 670       continue;
 671
 672     // kinematic cuts
 673     Double_t pt = mcP->Pt() ;
 674     if (pt<0.5)
 675       continue;
 676     Double_t eta = mcP->Eta();
 677     if (eta<-0.7||eta>0.7)
 678       continue;
 679     Double_t phi  = mcP->Phi();
 680     if (phi<1.0||phi>3.3)
 681       continue;
 682     // pion or eta meson or direct photon
 683     if(mcP->GetPdgCode() == 111) {
 684     } else if(mcP->GetPdgCode() == 221) {
 685     } else if(mcP->GetPdgCode() == 22 ) {
 686     } else
 687       continue;
 688     bool checkIfAlreadySaved = false;
 689     for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
 690       AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
 691       if(!mymc)
 692         continue;
 693       if(mymc->fLabel == iTrack)
 694         checkIfAlreadySaved = true;
 695     }
 696     if(!checkIfAlreadySaved)
 697       FillMcPart(mcP, ipart++, iTrack);
 698     for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
 699       if(id<=mcP->GetMother(0))
 700         continue;
 701       if(id<0 || id>nTracks)
 702         continue;
 703       TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
 704       if(!mcD)
 705         continue;
 706       FillMcPart(mcD, ipart++, id);
 707     }
 708   }
 709 }
 710
 711 //________________________________________________________________________
 712 void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
 713 {
 714   // Fill MC particles.
 715
 716   if(!mcP)
 717     return;
 718   TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
 719   AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
 720   mcp->fLabel = itrack ;
 721   mcp->fPdg = mcP->GetPdgCode() ;
 722   mcp->fPt = mcP->Pt() ;
 723   mcp->fEta = mcP->Eta() ;
 724   mcp->fPhi = mcP->Phi() ;
 725   mcp->fVR = vmcv.Perp();
 726   if(vmcv.Perp()>0){
 727     mcp->fVEta = vmcv.Eta() ;
 728     mcp->fVPhi = vmcv.Phi() ;
 729   }
 730   mcp->fMother = mcP->GetMother(0) ;
 731 }
 732
 733 //________________________________________________________________________
 734 Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
 735 {
 736   // Compute isolation based on tracks.
 737
 738   Double_t trkIsolation = 0;
 739   Double_t rad2 = radius*radius;
 740   Int_t ntrks = fSelPrimTracks->GetEntries();
 741   for(Int_t j = 0; j<ntrks; ++j) {
 742     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
 743     if (!track)
 744       continue;
 745     if (track->Pt()<pt)
 746       continue;
 747
 748     Float_t eta = track->Eta();
 749     Float_t phi = track->Phi();
 750     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
 751     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
 752     if(dist>rad2)
 753       continue;
 754     trkIsolation += track->Pt();
 755   }
 756   return trkIsolation;
 757 }
 758
 759 //________________________________________________________________________
 760 Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
 761 {
 762   // Get phi band.
 763
 764   if(!fSelPrimTracks)
 765     return 0;
 766   Int_t ntracks = fSelPrimTracks->GetEntries();
 767   Double_t loweta = eta - R;
 768   Double_t upeta = eta + R;
 769   Double_t upphi = phi + R;
 770   Double_t lowphi = phi - R;
 771   Double_t et = 0;
 772   for(Int_t itr=0; itr<ntracks; itr++){
 773     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
 774     if(!track)
 775       continue;
 776     if(track->Pt()<minpt)
 777       continue;
 778     if((track->Eta() < upeta) && (track->Eta() > loweta))
 779       continue;
 780     if((track->Phi() > upphi) || (track->Phi() < lowphi))
 781       continue;
 782     et+=track->Pt();
 783   }
 784   return et;
 785 }
 786
 787 //________________________________________________________________________
 788 Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 789 {
 790   // Calculate the energy of cross cells around the leading cell.
 791
 792   AliVCaloCells *cells = 0;
 793   cells = fESD->GetEMCALCells();
 794   if (!cells)
 795     return 0;
 796
 797   if (!fGeom)
 798     return 0;
 799
 800   Int_t iSupMod = -1;
 801   Int_t iTower  = -1;
 802   Int_t iIphi   = -1;
 803   Int_t iIeta   = -1;
 804   Int_t iphi    = -1;
 805   Int_t ieta    = -1;
 806   Int_t iphis   = -1;
 807   Int_t ietas   = -1;
 808
 809   Double_t crossEnergy = 0;
 810
 811   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 812   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 813
 814   Int_t ncells = cluster->GetNCells();
 815   for (Int_t i=0; i<ncells; i++) {
 816     Int_t cellAbsId = cluster->GetCellAbsId(i);
 817     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 818     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 819     Int_t aphidiff = TMath::Abs(iphi-iphis);
 820     if (aphidiff>1)
 821       continue;
 822     Int_t aetadiff = TMath::Abs(ieta-ietas);
 823     if (aetadiff>1)
 824       continue;
 825     if ( (aphidiff==1 && aetadiff==0) ||
 826         (aphidiff==0 && aetadiff==1) ) {
 827       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 828     }
 829   }
 830
 831   return crossEnergy;
 832 }
 833
 834 //________________________________________________________________________
 835 Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 836 {
 837   // Get maximum energy of attached cell.
 838
 839   id = -1;
 840
 841   AliVCaloCells *cells = 0;
 842   cells = fESD->GetEMCALCells();
 843   if (!cells)
 844     return 0;
 845
 846   Double_t maxe = 0;
 847   Int_t ncells = cluster->GetNCells();
 848   for (Int_t i=0; i<ncells; i++) {
 849     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 850     if (e>maxe) {
 851       maxe = e;
 852       id   = cluster->GetCellAbsId(i);
 853     }
 854   }
 855   return maxe;
 856 }
 857
 858 //________________________________________________________________________
 859 void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
 860 {
 861   // Called once at the end of the query
 862   if(fIsGrid)
 863     return;
 864   if (fTree) {
 865     TFile *f = OpenFile(1);
 866     TDirectory::TContext context(f);
 867     if (f)
 868       fTree->Write();
 869   }
 870 }