[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, 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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliAnalysisTaskSpectraAllChAOD class
//-----------------------------------------------------------------

#include "TChain.h"
#include "TTree.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TH2F.h"
#include "THnSparse.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliVParticle.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisTaskSpectraAllChAOD.h"
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisDataContainer.h"
#include "AliSpectraAODTrackCuts.h"
#include "AliSpectraAODEventCuts.h"
#include "AliHelperPID.h"
#include "AliPIDCombined.h"
#include "AliCentrality.h"
#include "TProof.h"
#include "AliVEvent.h"
#include "AliStack.h"
#include <TMCProcess.h>

#include <iostream>

using namespace AliHelperPIDNameSpace;
using namespace std;

ClassImp(AliAnalysisTaskSpectraAllChAOD)

//________________________________________________________________________
AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  
fAOD(0x0),
fTrackCuts(0x0),
fEventCuts(0x0),
fHelperPID(0x0),
fIsMC(0),
fDoDoubleCounting(0),
fFillOnlyEvents(0),
fCharge(0),
fVZEROside(0),
fOutput(0x0),
fnCentBins(20),
fnQvecBins(100),
fnNchBins(200),
fIsQvecCalibMode(0),
fQvecUpperLim(100),
fIsAOD160(1),
fnDCABins(60),
fDCAmin(-3),
fDCAmax(3),
fDCAzCut(999999.),
fQst(1),
fQtrk(0),
fQgenType(0),
fDoCentrSystCentrality(0)
{
  // Default constructor
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, AliSpectraAODEventCuts::Class());
  DefineOutput(3, AliSpectraAODTrackCuts::Class());
  DefineOutput(4, AliHelperPID::Class());
}

//________________________________________________________________________
void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
{
  // create output objects
  fOutput=new TList();
  fOutput->SetOwner();
  fOutput->SetName("fOutput");

  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
  if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");

  // binning common to all the THn
  const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
  const Int_t nptBins=34;
  //const Double_t ptBins[] = {0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
  //const Int_t nptBins=44;

  //dimensions of THnSparse for tracks
  const Int_t nvartrk=8;
  //                                             pt          cent          Q vec       IDrec          IDgen        isph         y             DCA 
  //  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2,        fnDCABins};
  //  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5,          fDCAmin};
  //  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5,         fDCAmax};
  //  THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
  //  NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
  //  NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
  //  NSparseHistTrk->SetBinEdges(0,ptBins);
  //  NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  //  NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  //  NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
  //  NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
  //  NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
  //  NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
  //  NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
  //  NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
  //  NSparseHistTrk->GetAxis(5)->SetTitle("isph");
  //  NSparseHistTrk->GetAxis(5)->SetName("isph");
  //  NSparseHistTrk->GetAxis(6)->SetTitle("y");
  //  NSparseHistTrk->GetAxis(6)->SetName("y");
  //  NSparseHistTrk->GetAxis(7)->SetTitle("dca");
  //  NSparseHistTrk->GetAxis(7)->SetName("dca");
  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2};
  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5};
  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5};
  THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
  NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
  NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
  NSparseHistTrk->SetBinEdges(0,ptBins);
  NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
  NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
  NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
  NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
  NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
  NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
  NSparseHistTrk->GetAxis(5)->SetTitle("isph");
  NSparseHistTrk->GetAxis(5)->SetName("isph");
  NSparseHistTrk->GetAxis(6)->SetTitle("y");
  NSparseHistTrk->GetAxis(6)->SetName("y");


  fOutput->Add(NSparseHistTrk);

  //dimensions of THnSparse for stack
  const Int_t nvarst=7;
  //                                             pt          cent    IDgen        isph        y    etaselected       Qvec gen
  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1,      fnQvecBins};
  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5,              0.};
  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5,      fQvecUpperLim};
  THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
  NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
  NSparseHistSt->SetBinEdges(0,ptBins);
  NSparseHistSt->GetAxis(0)->SetName("pT_rec");
  NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
  NSparseHistSt->GetAxis(2)->SetName("ID_gen");
  NSparseHistSt->GetAxis(3)->SetTitle("isph");
  NSparseHistSt->GetAxis(3)->SetName("isph");
  NSparseHistSt->GetAxis(4)->SetTitle("y");
  NSparseHistSt->GetAxis(4)->SetName("y");
  NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
  NSparseHistSt->GetAxis(5)->SetName("etaselected");
  NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
  NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
  fOutput->Add(NSparseHistSt);

  //dimensions of THnSparse for the normalization
  const Int_t nvarev=3;
  //                                             cent             Q vec                Nch
  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
  Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.};
  THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
  NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
  NSparseHistEv->GetAxis(1)->SetName("Q_vec");
  NSparseHistEv->GetAxis(2)->SetTitle("N charged");
  NSparseHistEv->GetAxis(2)->SetName("N_ch");
  fOutput->Add(NSparseHistEv);

  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//________________________________________________________________________

void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
{
  //Printf("An event");
  // main event loop
  fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
  if (!fAOD) {
    AliWarning("ERROR: AliAODEvent not available \n");
    return;
  }

  if (strcmp(fAOD->ClassName(), "AliAODEvent"))
  {
    AliFatal("Not processing AODs");
  }

  if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection

  //Default TPC priors
  if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME we should modify the task to change priors

  Double_t Qvec=0.;
  if(fIsQvecCalibMode){
    if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
    else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
    else if (fVZEROside==2)Qvec=fEventCuts->GetqTPC();
  }
  else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);

  Double_t QvecMC = 0.;
  if(fIsMC){
    if(fIsQvecCalibMode){
      QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside, fQgenType);
    }
    else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
  }

  Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();

  // First do MC to fill up the MC particle array
  TClonesArray *arrayMC = 0;
  if (fIsMC)
  {
    arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
    if (!arrayMC) {
      AliFatal("Error: MC particles branch not found!\n");
    }
    Int_t nMC = arrayMC->GetEntries();
    for (Int_t iMC = 0; iMC < nMC; iMC++)
    {
      AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
      if(!partMC->Charge()) continue;//Skip neutrals
      if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge

      //flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
      Double_t etaselected=-1.;
      if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;

      //pt     cent    IDgen        isph        y
      Double_t varSt[7];
      varSt[0]=partMC->Pt();
      varSt[1]=Cent;
      varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
      varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
      varSt[4]=partMC->Y();
      varSt[5]=etaselected;

      if(fQst==0) varSt[6]=Qvec;
      else varSt[6]=QvecMC;

      ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop

      //Printf("a particle");

    }
  }

  //main loop on tracks

  Int_t Nch = 0.;

  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
    AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
    if(!track) AliFatal("Not a standard AOD");
    if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge 
    if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
    if(!fFillOnlyEvents){
      Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
      Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
      Int_t IDgen=kSpUndefined;//set if MC
      Int_t isph=-999;
      //       Int_t iswd=-999;

      if (arrayMC) {
        AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
        if (!partMC) {
          AliError("Cannot get MC particle");
          continue;
        }
        IDgen=fHelperPID->GetParticleSpecies(partMC);
        isph=partMC->IsPhysicalPrimary();
        //iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014

        if(fIsAOD160){// enabled for new ADO160 only
          if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
          if(partMC->IsSecondaryFromMaterial()) isph=3.;

        }
      }

      /*** DCA ***/
      Double_t dcaxy = -999.;
      Double_t dcaz = -999.;

      Double_t p[2]; 
      if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}

      if(dcaz >= fDCAzCut) continue;

      //if the q vector is done using the TPC, we avoid overlap
      if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)continue;

      //pt     cent    Q vec     IDrec     IDgen       isph      y
      //      Double_t varTrk[8];
      Double_t varTrk[7];
      varTrk[0]=track->Pt();
      varTrk[1]=Cent;
      if(fIsMC && fQtrk==1) varTrk[2]=QvecMC;
      else varTrk[2]=Qvec;
      varTrk[3]=(Double_t)IDrec;
      varTrk[4]=(Double_t)IDgen;
      varTrk[5]=(Double_t)isph;
      varTrk[6]=y;
      //varTrk[7]=dcaxy;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop

      //for nsigma PID fill double counting of ID
      if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
        Bool_t *HasDC;
        HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
        for(Int_t ipart=0;ipart<kNSpecies;ipart++){
          if(HasDC[ipart]==kTRUE){
            varTrk[3]=(Double_t)ipart;
            ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
          }
        }
      }

      //fill all charged (3)
      varTrk[3]=3.;
      varTrk[4]=3.;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
    }//end if fFillOnlyEvents

    //Printf("a track");

    Nch++;
  } // end loop on tracks

  Double_t varEv[4];
  varEv[0]=Cent;
  varEv[1]=Qvec;
  varEv[2]=Nch;
  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop

  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//_________________________________________________________________
Bool_t  AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){

  //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
  //FIXME should update EventCuts?
  //FIXME add track->GetXYZ(p) method

  double xyz[3],cov[3];

  if (!trk->GetXYZ(xyz)) { // dca is not stored
    AliExternalTrackParam etp;
    etp.CopyFromVTrack(trk);
    AliVEvent* ev = (AliVEvent*)trk->GetEvent();
    if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
    if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
  }
  p[0] = xyz[0];
  p[1] = xyz[1];
  return kTRUE;

}

//_________________________________________________________________
void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
{
  // Terminate
}
Commit	Line	Data
c683985a	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, 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	// AliAnalysisTaskSpectraAllChAOD class
	18	//-----------------------------------------------------------------
	19
	20	#include "TChain.h"
	21	#include "TTree.h"
	22	#include "TLegend.h"
	23	#include "TH1F.h"
	24	#include "TH2F.h"
	25	#include "THnSparse.h"
	26	#include "TCanvas.h"
	27	#include "AliAnalysisTask.h"
	28	#include "AliAODTrack.h"
	29	#include "AliAODMCParticle.h"
	30	#include "AliVParticle.h"
	31	#include "AliAODEvent.h"
	32	#include "AliAODInputHandler.h"
	33	#include "AliAnalysisTaskSpectraAllChAOD.h"
	34	#include "AliAnalysisTaskESDfilter.h"
	35	#include "AliAnalysisDataContainer.h"
	36	#include "AliSpectraAODTrackCuts.h"
	37	#include "AliSpectraAODEventCuts.h"
	38	#include "AliHelperPID.h"
	39	#include "AliPIDCombined.h"
	40	#include "AliCentrality.h"
	41	#include "TProof.h"
	42	#include "AliVEvent.h"
	43	#include "AliStack.h"
	44	#include <TMCProcess.h>
	45
	46	#include <iostream>
	47
	48	using namespace AliHelperPIDNameSpace;
	49	using namespace std;
	50
	51	ClassImp(AliAnalysisTaskSpectraAllChAOD)
	52
	53	//________________________________________________________________________
	54	AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),
86bac2e4	55	fAOD(0x0),
	56	fTrackCuts(0x0),
	57	fEventCuts(0x0),
	58	fHelperPID(0x0),
	59	fIsMC(0),
	60	fDoDoubleCounting(0),
	61	fFillOnlyEvents(0),
	62	fCharge(0),
	63	fVZEROside(0),
	64	fOutput(0x0),
	65	fnCentBins(20),
	66	fnQvecBins(100),
	67	fnNchBins(200),
	68	fIsQvecCalibMode(0),
	69	fQvecUpperLim(100),
	70	fIsAOD160(1),
	71	fnDCABins(60),
	72	fDCAmin(-3),
	73	fDCAmax(3),
	74	fDCAzCut(999999.),
	75	fQst(1),
	76	fQtrk(0),
	77	fQgenType(0),
	78	fDoCentrSystCentrality(0)
c683985a	79	{
	80	// Default constructor
	81	DefineInput(0, TChain::Class());
	82	DefineOutput(1, TList::Class());
	83	DefineOutput(2, AliSpectraAODEventCuts::Class());
	84	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	85	DefineOutput(4, AliHelperPID::Class());
	86	}
	87
	88	//________________________________________________________________________
	89	void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
	90	{
	91	// create output objects
	92	fOutput=new TList();
	93	fOutput->SetOwner();
	94	fOutput->SetName("fOutput");
86bac2e4	95
c683985a	96	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
	97	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
	98	if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
86bac2e4	99
c683985a	100	// binning common to all the THn
c12a5f69	101	const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
	102	const Int_t nptBins=34;
	103	//const Double_t ptBins[] = {0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
	104	//const Int_t nptBins=44;
86bac2e4	105
c683985a	106	//dimensions of THnSparse for tracks
3dd5ed48	107	const Int_t nvartrk=8;
1fdecdaf	108	// pt cent Q vec IDrec IDgen isph y DCA
86bac2e4	109	// Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 3, 2, fnDCABins};
	110	// Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, 0.5, -0.5, fDCAmin};
	111	// Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 3.5, 0.5, fDCAmax};
	112	// THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
	113	// NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
	114	// NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
	115	// NSparseHistTrk->SetBinEdges(0,ptBins);
	116	// NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	117	// NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	118	// NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
	119	// NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
	120	// NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
	121	// NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
	122	// NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
	123	// NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
	124	// NSparseHistTrk->GetAxis(5)->SetTitle("isph");
	125	// NSparseHistTrk->GetAxis(5)->SetName("isph");
	126	// NSparseHistTrk->GetAxis(6)->SetTitle("y");
	127	// NSparseHistTrk->GetAxis(6)->SetName("y");
	128	// NSparseHistTrk->GetAxis(7)->SetTitle("dca");
	129	// NSparseHistTrk->GetAxis(7)->SetName("dca");
	130	Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 3, 2};
	131	Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, 0.5, -0.5};
	132	Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 3.5, 0.5};
c683985a	133	THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
	134	NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
	135	NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
	136	NSparseHistTrk->SetBinEdges(0,ptBins);
	137	NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	138	NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	139	NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
	140	NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
	141	NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
	142	NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
	143	NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
	144	NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
	145	NSparseHistTrk->GetAxis(5)->SetTitle("isph");
	146	NSparseHistTrk->GetAxis(5)->SetName("isph");
eba59d12	147	NSparseHistTrk->GetAxis(6)->SetTitle("y");
eba59d12	148	NSparseHistTrk->GetAxis(6)->SetName("y");
86bac2e4	149
86bac2e4	150
c683985a	151	fOutput->Add(NSparseHistTrk);
86bac2e4	152
c683985a	153	//dimensions of THnSparse for stack
ed33b9ef	154	const Int_t nvarst=7;
	155	// pt cent IDgen isph y etaselected Qvec gen
	156	Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2, 1, fnQvecBins};
	157	Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5, 0.5, 0.};
	158	Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5, 1.5, fQvecUpperLim};
c683985a	159	THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
	160	NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
	161	NSparseHistSt->SetBinEdges(0,ptBins);
	162	NSparseHistSt->GetAxis(0)->SetName("pT_rec");
	163	NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	164	NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	165	NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
	166	NSparseHistSt->GetAxis(2)->SetName("ID_gen");
	167	NSparseHistSt->GetAxis(3)->SetTitle("isph");
	168	NSparseHistSt->GetAxis(3)->SetName("isph");
	169	NSparseHistSt->GetAxis(4)->SetTitle("y");
	170	NSparseHistSt->GetAxis(4)->SetName("y");
eba59d12	171	NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
eba59d12	172	NSparseHistSt->GetAxis(5)->SetName("etaselected");
ed33b9ef	173	NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
ed33b9ef	174	NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
c683985a	175	fOutput->Add(NSparseHistSt);
86bac2e4	176
c683985a	177	//dimensions of THnSparse for the normalization
86bac2e4	178	const Int_t nvarev=3;
	179	// cent Q vec Nch
	180	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins};
	181	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
	182	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 2000.};
c683985a	183	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	184	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	185	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	186	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	187	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
e19bc133	188	NSparseHistEv->GetAxis(2)->SetTitle("N charged");
e19bc133	189	NSparseHistEv->GetAxis(2)->SetName("N_ch");
c683985a	190	fOutput->Add(NSparseHistEv);
86bac2e4	191
c683985a	192	PostData(1, fOutput );
	193	PostData(2, fEventCuts);
	194	PostData(3, fTrackCuts);
	195	PostData(4, fHelperPID);
	196	}
	197
	198	//________________________________________________________________________
	199
	200	void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
	201	{
	202	//Printf("An event");
	203	// main event loop
	204	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	205	if (!fAOD) {
	206	AliWarning("ERROR: AliAODEvent not available \n");
	207	return;
	208	}
86bac2e4	209
c683985a	210	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
86bac2e4	211	{
	212	AliFatal("Not processing AODs");
	213	}
	214
c683985a	215	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
	216
	217	//Default TPC priors
13dab2da	218	if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME we should modify the task to change priors
86bac2e4	219
ed33b9ef	220	Double_t Qvec=0.;
	221	if(fIsQvecCalibMode){
	222	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	223	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
a72cae45	224	else if (fVZEROside==2)Qvec=fEventCuts->GetqTPC();
ed33b9ef	225	}
ed33b9ef	226	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
86bac2e4	227
ed33b9ef	228	Double_t QvecMC = 0.;
ed33b9ef	229	if(fIsMC){
eba59d12	230	if(fIsQvecCalibMode){
725da720	231	QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside, fQgenType);
eba59d12	232	}
725da720	233	else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
c683985a	234	}
86bac2e4	235
509a25f9	236	Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
86bac2e4	237
c683985a	238	// First do MC to fill up the MC particle array
	239	TClonesArray *arrayMC = 0;
	240	if (fIsMC)
86bac2e4	241	{
	242	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	243	if (!arrayMC) {
	244	AliFatal("Error: MC particles branch not found!\n");
	245	}
	246	Int_t nMC = arrayMC->GetEntries();
	247	for (Int_t iMC = 0; iMC < nMC; iMC++)
c683985a	248	{
86bac2e4	249	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);
	250	if(!partMC->Charge()) continue;//Skip neutrals
	251	if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
	252
	253	//flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
	254	Double_t etaselected=-1.;
	255	if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
	256
	257	//pt cent IDgen isph y
	258	Double_t varSt[7];
	259	varSt[0]=partMC->Pt();
	260	varSt[1]=Cent;
	261	varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	262	varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	263	varSt[4]=partMC->Y();
	264	varSt[5]=etaselected;
	265
	266	if(fQst==0) varSt[6]=Qvec;
	267	else varSt[6]=QvecMC;
	268
	269	((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
	270
	271	//Printf("a particle");
	272
c683985a	273	}
86bac2e4	274	}
86bac2e4	275
c683985a	276	//main loop on tracks
86bac2e4	277
e19bc133	278	Int_t Nch = 0.;
86bac2e4	279
c683985a	280	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
f15c1f69	281	AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
f15c1f69	282	if(!track) AliFatal("Not a standard AOD");
c683985a	283	if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
c683985a	284	if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
dbbc07ac	285	if(!fFillOnlyEvents){
	286	Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
	287	Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
	288	Int_t IDgen=kSpUndefined;//set if MC
	289	Int_t isph=-999;
86bac2e4	290	// Int_t iswd=-999;
86bac2e4	291
dbbc07ac	292	if (arrayMC) {
86bac2e4	293	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));
	294	if (!partMC) {
	295	AliError("Cannot get MC particle");
	296	continue;
	297	}
	298	IDgen=fHelperPID->GetParticleSpecies(partMC);
	299	isph=partMC->IsPhysicalPrimary();
	300	//iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
	301
	302	if(fIsAOD160){// enabled for new ADO160 only
	303	if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
	304	if(partMC->IsSecondaryFromMaterial()) isph=3.;
	305
	306	}
c683985a	307	}
86bac2e4	308
93c24428	309	/* DCA */
93c24428	310	Double_t dcaxy = -999.;
3dd5ed48	311	Double_t dcaz = -999.;
86bac2e4	312
665f3db3	313	Double_t p[2];
3dd5ed48	314	if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
86bac2e4	315
3dd5ed48	316	if(dcaz >= fDCAzCut) continue;
896040ab	317
	318	//if the q vector is done using the TPC, we avoid overlap
	319	if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)continue;
86bac2e4	320
	321	//pt cent Q vec IDrec IDgen isph y
	322	// Double_t varTrk[8];
	323	Double_t varTrk[7];
dbbc07ac	324	varTrk[0]=track->Pt();
dbbc07ac	325	varTrk[1]=Cent;
09f5725e	326	if(fIsMC && fQtrk==1) varTrk[2]=QvecMC;
86bac2e4	327	else varTrk[2]=Qvec;
dbbc07ac	328	varTrk[3]=(Double_t)IDrec;
	329	varTrk[4]=(Double_t)IDgen;
	330	varTrk[5]=(Double_t)isph;
eba59d12	331	varTrk[6]=y;
86bac2e4	332	//varTrk[7]=dcaxy;
dbbc07ac	333	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
86bac2e4	334
dbbc07ac	335	//for nsigma PID fill double counting of ID
dbbc07ac	336	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
86bac2e4	337	Bool_t *HasDC;
	338	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	339	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	340	if(HasDC[ipart]==kTRUE){
	341	varTrk[3]=(Double_t)ipart;
	342	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	343	}
	344	}
c683985a	345	}
86bac2e4	346
dbbc07ac	347	//fill all charged (3)
	348	varTrk[3]=3.;
	349	varTrk[4]=3.;
	350	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	351	}//end if fFillOnlyEvents
86bac2e4	352
c683985a	353	//Printf("a track");
86bac2e4	354
dbbc07ac	355	Nch++;
c683985a	356	} // end loop on tracks
86bac2e4	357
ed33b9ef	358	Double_t varEv[4];
c683985a	359	varEv[0]=Cent;
c683985a	360	varEv[1]=Qvec;
e19bc133	361	varEv[2]=Nch;
c683985a	362	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
86bac2e4	363
c683985a	364	PostData(1, fOutput );
	365	PostData(2, fEventCuts);
	366	PostData(3, fTrackCuts);
	367	PostData(4, fHelperPID);
	368	}
	369
665f3db3	370	//_________________________________________________________________
665f3db3	371	Bool_t AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
86bac2e4	372
665f3db3	373	//AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
	374	//FIXME should update EventCuts?
	375	//FIXME add track->GetXYZ(p) method
86bac2e4	376
665f3db3	377	double xyz[3],cov[3];
86bac2e4	378
665f3db3	379	if (!trk->GetXYZ(xyz)) { // dca is not stored
	380	AliExternalTrackParam etp;
	381	etp.CopyFromVTrack(trk);
	382	AliVEvent* ev = (AliVEvent*)trk->GetEvent();
	383	if (!ev) {/printf("Event is not connected to the track\n");/ return kFALSE;}
	384	if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
	385	}
	386	p[0] = xyz[0];
	387	p[1] = xyz[1];
	388	return kTRUE;
	389
	390	}
	391
c683985a	392	//_________________________________________________________________
	393	void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
	394	{
	395	// Terminate
	396	}