[u/mrichter/AliRoot.git] / PWGJE / PWGJE / AliPWG4CosmicCandidates.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------
// Analysis task looking for cosmic candidates embedded in a p-p event
// Task checks if particles are back-to-back in eta and phi
//             
//
// Author : Marta Verweij - UU - marta.verweij@cern.ch
//-----------------------------------------------------------------------

#include "TVector3.h"
#include <iostream>
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TList.h"
#include "TChain.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"

#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliExternalTrackParam.h"

#include "AliLog.h"

#include "AliPWG4CosmicCandidates.h"

//using namespace std; //required for resolving the 'cout' symbol
using namespace std;

ClassImp(AliPWG4CosmicCandidates)

//________________________________________________________________________
AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
: AliAnalysisTaskSE(),
  fTrackCuts(0), 
  fPtMin(5.),
  fMaxCosmicAngle(0.002),
  fNEventAll(0),
  fNEventSel(0),
  fPtSignedCosmicCandidates(0),
  fDeltaPtCosmicCandidates(0),
  fDeltaPhiSumEta(0),
  fDCAZCosmicCandidates(0),
  fDCARCosmicCandidates(0),
  fTheta(0),
  fThetaZoom(0),
  fThetaPt1Pt2(0),
  fThetaPt1Pt2Signed(0),
  fDeltaPhiSumEtaPt1(0),
  fDeltaPhiSumEtaPt2(0),
  fThetaDCAZ1DCAZ2(0),
  fRisol(0),
  fRisolTheta(0),
  fHistListCosmics(0)
{
  //
  // Default constructor
  //
}

//________________________________________________________________________
AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
  : AliAnalysisTaskSE(name),
    fTrackCuts(0), 
    fPtMin(5.),
    fMaxCosmicAngle(0.002),
    fNEventAll(0),
    fNEventSel(0),
    fPtSignedCosmicCandidates(0),
    fDeltaPtCosmicCandidates(0),
    fDeltaPhiSumEta(0),
    fDCAZCosmicCandidates(0),
    fDCARCosmicCandidates(0),
    fTheta(0),
    fThetaZoom(0),
    fThetaPt1Pt2(0),
    fThetaPt1Pt2Signed(0),
    fDeltaPhiSumEtaPt1(0),
    fDeltaPhiSumEtaPt2(0),
    fThetaDCAZ1DCAZ2(0),
    fRisol(0),
    fRisolTheta(0),
    fHistListCosmics(0)
{
  // Constructor. Initialization of Inputs and Outputs

  // Define input and output slots here
  // Input slot #0 works with a TChain
  //  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TList
  DefineOutput(1, TList::Class());
  // Output slot #2 writes into a AliESDtrackCuts
  DefineOutput(2, AliESDtrackCuts::Class());
}
//________________________________________________________________________
AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
  : AliAnalysisTaskSE(res),
    fTrackCuts(0), 
    fPtMin(5.),
    fMaxCosmicAngle(0.002),
    fNEventAll(0),
    fNEventSel(0),
    fPtSignedCosmicCandidates(0),
    fDeltaPtCosmicCandidates(0),
    fDeltaPhiSumEta(0),
    fDCAZCosmicCandidates(0),
    fDCARCosmicCandidates(0),
    fTheta(0),
    fThetaZoom(0),
    fThetaPt1Pt2(0),
    fThetaPt1Pt2Signed(0),
    fDeltaPhiSumEtaPt1(0),
    fDeltaPhiSumEtaPt2(0),
    fThetaDCAZ1DCAZ2(0),
    fRisol(0),
    fRisolTheta(0),
    fHistListCosmics(0)
{
  //
  // Dummy copy constructor
  //
}

//________________________________________________________________________
AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
{
  //
  // Dummy assignment operator
  //
  return *this;
}

//________________________________________________________________________
void AliPWG4CosmicCandidates::LocalInit()
{
  //
  // Only called once at beginning
  //
  PostData(2,fTrackCuts);
}

//________________________________________________________________________
void AliPWG4CosmicCandidates::UserCreateOutputObjects()
{
  //
  // Create output objects
  // Called once
  AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n")); 

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE); 
  
  OpenFile(1);
  fHistListCosmics = new TList();

  fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
  fHistListCosmics->Add(fNEventAll);
  fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
  fHistListCosmics->Add(fNEventSel);

  Float_t fgkPtMin=0.;
  Float_t fgkPtMax=100.;
  Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
  Int_t fgkNPtBins2D= (int)((fgkPtMax-fgkPtMin)/4.);

  Int_t fgkNPhiBins=18;
  Float_t kMinPhi = -0.5*TMath::Pi();
  Float_t kMaxPhi = 3./2.*TMath::Pi();

  Int_t fgkNThetaBins=fgkNPhiBins*8;
  Float_t kMinTheta = -0.5*TMath::Pi();
  Float_t kMaxTheta = 3./2.*TMath::Pi();

  Int_t fgkNDCARBins=40;
  Float_t fgkDCARMin = -0.2;
  Float_t fgkDCARMax = 0.2;
  Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
  for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;

  Int_t fgkNDCAZBins=40;
  Float_t fgkDCAZMin = -2.;
  Float_t fgkDCAZMax = 2.;
  Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
  for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;

  fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*fgkNPtBins, -1.*fgkPtMax, fgkPtMax);
  fHistListCosmics->Add(fPtSignedCosmicCandidates);  

  fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
  fHistListCosmics->Add(fDeltaPtCosmicCandidates);  

  fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
  fHistListCosmics->Add(fDeltaPhiSumEta);  

  fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
  fHistListCosmics->Add(fDCAZCosmicCandidates);

  fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
  fHistListCosmics->Add(fDCARCosmicCandidates);

  fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
  fHistListCosmics->Add(fTheta);

  fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
  fHistListCosmics->Add(fThetaZoom);

  fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,fgkPtMin,fgkPtMax,fgkNPtBins2D,fgkPtMin,fgkPtMax);
  fHistListCosmics->Add(fThetaPt1Pt2);

  fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax);
  fHistListCosmics->Add(fThetaPt1Pt2Signed);

  fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
  fHistListCosmics->Add(fDeltaPhiSumEtaPt1);

  fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
  fHistListCosmics->Add(fDeltaPhiSumEtaPt2);

  fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
  fHistListCosmics->Add(fThetaDCAZ1DCAZ2);

  fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
  fHistListCosmics->Add(fRisol);

  fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
  fHistListCosmics->Add(fRisolTheta);

  TH1::AddDirectory(oldStatus); 

  PostData(1, fHistListCosmics);  

  if(binsDCAR) delete [] binsDCAR;
  if(binsDCAZ) delete [] binsDCAZ;

}

//________________________________________________________________________
void AliPWG4CosmicCandidates::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  //

  // All events without selection
  fNEventAll->Fill(0.);
   
  if (!fInputEvent) {
    AliDebug(2,Form("ERROR: fESD not available"));
    cout << "ERROR: fESD not available" << endl;
    PostData(1, fHistListCosmics);
    return;
  }

  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
  if(!vtx){
    PostData(1, fHistListCosmics);
    return;
  }

// Need vertex cut
  TString vtxName(vtx->GetName());
  if( vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
    // Post output data
    PostData(1, fHistListCosmics);
    return;
  }

  //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
  double primVtx[3];
  vtx->GetXYZ(primVtx);
  if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
    // Post output data
    PostData(1, fHistListCosmics);
    return;
  }
  if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){ 
    // Post output data
    PostData(1, fHistListCosmics);
    return;
  }
  Int_t nTracks = fInputEvent->GetNumberOfTracks();

  if(!fTrackCuts) {
   // Post output data
    PostData(1, fHistListCosmics);
    return;
  }

  fNEventSel->Fill(0.);

  Float_t dcaR[2] = {0.,0.};
  Float_t dcaZ[2] = {0.,0.};

  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {

    AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
    if (!track1)  continue;
    if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
    if(track1->Pt()<fPtMin) continue;
    //Start 2nd track loop to look for correlations
    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
      AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
      if(!track2) continue;
      if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
         
      //Check if back-to-back
      Double_t mom1[3],mom2[3];
      track1->GetPxPyPz(mom1);
      track2->GetPxPyPz(mom2);
      //     Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
      TVector3 momv1(mom1[0],mom1[1],mom1[2]);
      TVector3 momv2(mom2[0],mom2[1],mom2[2]);
      //Double_t theta = momv1.Angle(momv2);
      Double_t theta = momv1.Phi()-momv2.Phi();
      if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
    
      fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
      Float_t deltaPhi = track1->Phi()-track2->Phi();
      if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
      fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());

      track1->GetImpactParameters(dcaR[0],dcaZ[0]);
      track2->GetImpactParameters(dcaR[1],dcaZ[1]);

      if(track2->Pt()<0.5) continue;
      Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
      fRisol->Fill(rIsol); //Fill R histogram
      if(track2->Pt()<fPtMin) continue;

	fTheta->Fill(theta);
	fThetaZoom->Fill(theta);
	fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
	fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
	fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
	fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
	fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
	fRisolTheta->Fill(rIsol,theta);
	if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
	  fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
	  fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
	  fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
	}

    } // track2 loop
    
  } // track1 loop 

  PostData(1,fHistListCosmics);
}      

//________________________________________________________________________
void AliPWG4CosmicCandidates::Terminate(Option_t *) 
{
  //
  // Called once at the end of the query
  //
  
}
Commit	Line	Data
	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------------
	17	// Analysis task looking for cosmic candidates embedded in a p-p event
	18	// Task checks if particles are back-to-back in eta and phi
	19	//
	20	//
	21	// Author : Marta Verweij - UU - marta.verweij@cern.ch
	22	//-----------------------------------------------------------------------
	23
	24	#include "TVector3.h"
	25	#include <iostream>
	26	#include "TH1.h"
	27	#include "TH2.h"
	28	#include "TH3.h"
	29	#include "TList.h"
	30	#include "TChain.h"
	31
	32	#include "AliAnalysisTask.h"
	33	#include "AliAnalysisManager.h"
	34
	35	#include "AliESDEvent.h"
	36	#include "AliESDInputHandler.h"
	37
	38	#include "AliESDtrack.h"
	39	#include "AliESDtrackCuts.h"
	40	#include "AliExternalTrackParam.h"
	41
	42	#include "AliLog.h"
	43
	44	#include "AliPWG4CosmicCandidates.h"
	45
	46	//using namespace std; //required for resolving the 'cout' symbol
	47	using namespace std;
	48
	49	ClassImp(AliPWG4CosmicCandidates)
	50
	51	//________________________________________________________________________
	52	AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
	53	: AliAnalysisTaskSE(),
	54	fTrackCuts(0),
	55	fPtMin(5.),
	56	fMaxCosmicAngle(0.002),
	57	fNEventAll(0),
	58	fNEventSel(0),
	59	fPtSignedCosmicCandidates(0),
	60	fDeltaPtCosmicCandidates(0),
	61	fDeltaPhiSumEta(0),
	62	fDCAZCosmicCandidates(0),
	63	fDCARCosmicCandidates(0),
	64	fTheta(0),
	65	fThetaZoom(0),
	66	fThetaPt1Pt2(0),
	67	fThetaPt1Pt2Signed(0),
	68	fDeltaPhiSumEtaPt1(0),
	69	fDeltaPhiSumEtaPt2(0),
	70	fThetaDCAZ1DCAZ2(0),
	71	fRisol(0),
	72	fRisolTheta(0),
	73	fHistListCosmics(0)
	74	{
	75	//
	76	// Default constructor
	77	//
	78	}
	79
	80	//________________________________________________________________________
	81	AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
	82	: AliAnalysisTaskSE(name),
	83	fTrackCuts(0),
	84	fPtMin(5.),
	85	fMaxCosmicAngle(0.002),
	86	fNEventAll(0),
	87	fNEventSel(0),
	88	fPtSignedCosmicCandidates(0),
	89	fDeltaPtCosmicCandidates(0),
	90	fDeltaPhiSumEta(0),
	91	fDCAZCosmicCandidates(0),
	92	fDCARCosmicCandidates(0),
	93	fTheta(0),
	94	fThetaZoom(0),
	95	fThetaPt1Pt2(0),
	96	fThetaPt1Pt2Signed(0),
	97	fDeltaPhiSumEtaPt1(0),
	98	fDeltaPhiSumEtaPt2(0),
	99	fThetaDCAZ1DCAZ2(0),
	100	fRisol(0),
	101	fRisolTheta(0),
	102	fHistListCosmics(0)
	103	{
	104	// Constructor. Initialization of Inputs and Outputs
	105
	106	// Define input and output slots here
	107	// Input slot #0 works with a TChain
	108	// DefineInput(0, TChain::Class());
	109	// Output slot #0 id reserved by the base class for AOD
	110	// Output slot #1 writes into a TList
	111	DefineOutput(1, TList::Class());
	112	// Output slot #2 writes into a AliESDtrackCuts
	113	DefineOutput(2, AliESDtrackCuts::Class());
	114	}
	115	//________________________________________________________________________
	116	AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
	117	: AliAnalysisTaskSE(res),
	118	fTrackCuts(0),
	119	fPtMin(5.),
	120	fMaxCosmicAngle(0.002),
	121	fNEventAll(0),
	122	fNEventSel(0),
	123	fPtSignedCosmicCandidates(0),
	124	fDeltaPtCosmicCandidates(0),
	125	fDeltaPhiSumEta(0),
	126	fDCAZCosmicCandidates(0),
	127	fDCARCosmicCandidates(0),
	128	fTheta(0),
	129	fThetaZoom(0),
	130	fThetaPt1Pt2(0),
	131	fThetaPt1Pt2Signed(0),
	132	fDeltaPhiSumEtaPt1(0),
	133	fDeltaPhiSumEtaPt2(0),
	134	fThetaDCAZ1DCAZ2(0),
	135	fRisol(0),
	136	fRisolTheta(0),
	137	fHistListCosmics(0)
	138	{
	139	//
	140	// Dummy copy constructor
	141	//
	142	}
	143
	144	//________________________________________________________________________
	145	AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /trclass/)
	146	{
	147	//
	148	// Dummy assignment operator
	149	//
	150	return *this;
	151	}
	152
	153	//________________________________________________________________________
	154	void AliPWG4CosmicCandidates::LocalInit()
	155	{
	156	//
	157	// Only called once at beginning
	158	//
	159	PostData(2,fTrackCuts);
	160	}
	161
	162	//________________________________________________________________________
	163	void AliPWG4CosmicCandidates::UserCreateOutputObjects()
	164	{
	165	//
	166	// Create output objects
	167	// Called once
	168	AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n"));
	169
	170	Bool_t oldStatus = TH1::AddDirectoryStatus();
	171	TH1::AddDirectory(kFALSE);
	172
	173	OpenFile(1);
	174	fHistListCosmics = new TList();
	175
	176	fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
	177	fHistListCosmics->Add(fNEventAll);
	178	fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
	179	fHistListCosmics->Add(fNEventSel);
	180
	181	Float_t fgkPtMin=0.;
	182	Float_t fgkPtMax=100.;
	183	Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
	184	Int_t fgkNPtBins2D= (int)((fgkPtMax-fgkPtMin)/4.);
	185
	186	Int_t fgkNPhiBins=18;
	187	Float_t kMinPhi = -0.5*TMath::Pi();
	188	Float_t kMaxPhi = 3./2.*TMath::Pi();
	189
	190	Int_t fgkNThetaBins=fgkNPhiBins*8;
	191	Float_t kMinTheta = -0.5*TMath::Pi();
	192	Float_t kMaxTheta = 3./2.*TMath::Pi();
	193
	194	Int_t fgkNDCARBins=40;
	195	Float_t fgkDCARMin = -0.2;
	196	Float_t fgkDCARMax = 0.2;
	197	Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
	198	for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
	199
	200	Int_t fgkNDCAZBins=40;
	201	Float_t fgkDCAZMin = -2.;
	202	Float_t fgkDCAZMax = 2.;
	203	Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
	204	for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
	205
	206	fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2fgkNPtBins, -1.fgkPtMax, fgkPtMax);
	207	fHistListCosmics->Add(fPtSignedCosmicCandidates);
	208
	209	fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
	210	fHistListCosmics->Add(fDeltaPtCosmicCandidates);
	211
	212	fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
	213	fHistListCosmics->Add(fDeltaPhiSumEta);
	214
	215	fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
	216	fHistListCosmics->Add(fDCAZCosmicCandidates);
	217
	218	fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
	219	fHistListCosmics->Add(fDCARCosmicCandidates);
	220
	221	fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
	222	fHistListCosmics->Add(fTheta);
	223
	224	fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
	225	fHistListCosmics->Add(fThetaZoom);
	226
	227	fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,fgkPtMin,fgkPtMax,fgkNPtBins2D,fgkPtMin,fgkPtMax);
	228	fHistListCosmics->Add(fThetaPt1Pt2);
	229
	230	fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,-1.fgkPtMax,fgkPtMax,fgkNPtBins2D,-1.fgkPtMax,fgkPtMax);
	231	fHistListCosmics->Add(fThetaPt1Pt2Signed);
	232
	233	fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
	234	fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
	235
	236	fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
	237	fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
	238
	239	fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
	240	fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
	241
	242	fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
	243	fHistListCosmics->Add(fRisol);
	244
	245	fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
	246	fHistListCosmics->Add(fRisolTheta);
	247
	248	TH1::AddDirectory(oldStatus);
	249
	250	PostData(1, fHistListCosmics);
	251
	252	if(binsDCAR) delete [] binsDCAR;
	253	if(binsDCAZ) delete [] binsDCAZ;
	254
	255	}
	256
	257	//________________________________________________________________________
	258	void AliPWG4CosmicCandidates::UserExec(Option_t *)
	259	{
	260	// Main loop
	261	// Called for each event
	262	//
	263
	264	// All events without selection
	265	fNEventAll->Fill(0.);
	266
	267	if (!fInputEvent) {
	268	AliDebug(2,Form("ERROR: fESD not available"));
	269	cout << "ERROR: fESD not available" << endl;
	270	PostData(1, fHistListCosmics);
	271	return;
	272	}
	273
	274	const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
	275	if(!vtx){
	276	PostData(1, fHistListCosmics);
	277	return;
	278	}
	279
	280	// Need vertex cut
	281	TString vtxName(vtx->GetName());
	282	if( vtx->GetNContributors() < 2 \|\| (vtxName.Contains("TPCVertex")) ) {
	283	// Post output data
	284	PostData(1, fHistListCosmics);
	285	return;
	286	}
	287
	288	// AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
	289	double primVtx[3];
	290	vtx->GetXYZ(primVtx);
	291	if(TMath::Sqrt(primVtx[0]primVtx[0] + primVtx[1]primVtx[1])>1. \|\| TMath::Abs(primVtx[2]>10.)){
	292	// Post output data
	293	PostData(1, fHistListCosmics);
	294	return;
	295	}
	296	if(!fInputEvent->GetNumberOfTracks() \|\| fInputEvent->GetNumberOfTracks()<2){
	297	// Post output data
	298	PostData(1, fHistListCosmics);
	299	return;
	300	}
	301	Int_t nTracks = fInputEvent->GetNumberOfTracks();
	302
	303	if(!fTrackCuts) {
	304	// Post output data
	305	PostData(1, fHistListCosmics);
	306	return;
	307	}
	308
	309	fNEventSel->Fill(0.);
	310
	311	Float_t dcaR[2] = {0.,0.};
	312	Float_t dcaZ[2] = {0.,0.};
	313
	314	for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
	315
	316	AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
	317	if (!track1) continue;
	318	if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
	319	if(track1->Pt()<fPtMin) continue;
	320	//Start 2nd track loop to look for correlations
	321	for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
	322	AliESDtrack track2 = (AliESDtrack)fInputEvent->GetTrack(iTrack2);
	323	if(!track2) continue;
	324	if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
	325
	326	//Check if back-to-back
	327	Double_t mom1[3],mom2[3];
	328	track1->GetPxPyPz(mom1);
	329	track2->GetPxPyPz(mom2);
	330	// Double_t cosTheta = (mom1[0]mom2[0]+mom1[1]mom2[1]+mom1[2]mom2[2])/( TMath::Sqrt(mom1[0]mom1[0]+mom1[1]mom1[1]+mom1[2]mom1[2])TMath::Sqrt(mom2[0]mom2[0]+mom2[1]mom2[1]+mom2[2]mom2[2]) );
	331	TVector3 momv1(mom1[0],mom1[1],mom1[2]);
	332	TVector3 momv2(mom2[0],mom2[1],mom2[2]);
	333	//Double_t theta = momv1.Angle(momv2);
	334	Double_t theta = momv1.Phi()-momv2.Phi();
	335	if(theta<-0.5TMath::Pi()) theta+=2.TMath::Pi();
	336
	337	fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
	338	Float_t deltaPhi = track1->Phi()-track2->Phi();
	339	if(deltaPhi<-0.5TMath::Pi()) deltaPhi+=2.TMath::Pi();
	340	fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
	341
	342	track1->GetImpactParameters(dcaR[0],dcaZ[0]);
	343	track2->GetImpactParameters(dcaR[1],dcaZ[1]);
	344
	345	if(track2->Pt()<0.5) continue;
	346	Double_t rIsol = TMath::Sqrt( deltaPhideltaPhi+(track1->Eta()-track2->Eta())(track1->Eta()-track2->Eta()) );
	347	fRisol->Fill(rIsol); //Fill R histogram
	348	if(track2->Pt()<fPtMin) continue;
	349
	350	fTheta->Fill(theta);
	351	fThetaZoom->Fill(theta);
	352	fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
	353	fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()track1->Pt(),track2->GetSign()track2->Pt());
	354	fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
	355	fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
	356	fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
	357	fRisolTheta->Fill(rIsol,theta);
	358	if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
	359	fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
	360	fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
	361	fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
	362	}
	363
	364	} // track2 loop
	365
	366	} // track1 loop
	367
	368	PostData(1,fHistListCosmics);
	369	}
	370
	371	//________________________________________________________________________
	372	void AliPWG4CosmicCandidates::Terminate(Option_t *)
	373	{
	374	//
	375	// Called once at the end of the query
	376	//
	377
	378	}