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