[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");
  const Int_t nvartrk=7;
  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
494aafcd	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");
494aafcd	130	const Int_t nvartrk=7;
86bac2e4	131	Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 3, 2};
	132	Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, 0.5, -0.5};
	133	Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 3.5, 0.5};
c683985a	134	THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
	135	NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
	136	NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
	137	NSparseHistTrk->SetBinEdges(0,ptBins);
	138	NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	139	NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	140	NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
	141	NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
	142	NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
	143	NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
	144	NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
	145	NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
	146	NSparseHistTrk->GetAxis(5)->SetTitle("isph");
	147	NSparseHistTrk->GetAxis(5)->SetName("isph");
eba59d12	148	NSparseHistTrk->GetAxis(6)->SetTitle("y");
eba59d12	149	NSparseHistTrk->GetAxis(6)->SetName("y");
86bac2e4	150
86bac2e4	151
c683985a	152	fOutput->Add(NSparseHistTrk);
86bac2e4	153
c683985a	154	//dimensions of THnSparse for stack
ed33b9ef	155	const Int_t nvarst=7;
	156	// pt cent IDgen isph y etaselected Qvec gen
	157	Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2, 1, fnQvecBins};
	158	Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5, 0.5, 0.};
	159	Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5, 1.5, fQvecUpperLim};
c683985a	160	THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
	161	NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
	162	NSparseHistSt->SetBinEdges(0,ptBins);
	163	NSparseHistSt->GetAxis(0)->SetName("pT_rec");
	164	NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	165	NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	166	NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
	167	NSparseHistSt->GetAxis(2)->SetName("ID_gen");
	168	NSparseHistSt->GetAxis(3)->SetTitle("isph");
	169	NSparseHistSt->GetAxis(3)->SetName("isph");
	170	NSparseHistSt->GetAxis(4)->SetTitle("y");
	171	NSparseHistSt->GetAxis(4)->SetName("y");
eba59d12	172	NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
eba59d12	173	NSparseHistSt->GetAxis(5)->SetName("etaselected");
ed33b9ef	174	NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
ed33b9ef	175	NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
c683985a	176	fOutput->Add(NSparseHistSt);
86bac2e4	177
c683985a	178	//dimensions of THnSparse for the normalization
86bac2e4	179	const Int_t nvarev=3;
	180	// cent Q vec Nch
	181	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins};
	182	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
	183	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 2000.};
c683985a	184	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	185	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	186	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	187	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	188	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
e19bc133	189	NSparseHistEv->GetAxis(2)->SetTitle("N charged");
e19bc133	190	NSparseHistEv->GetAxis(2)->SetName("N_ch");
c683985a	191	fOutput->Add(NSparseHistEv);
86bac2e4	192
c683985a	193	PostData(1, fOutput );
	194	PostData(2, fEventCuts);
	195	PostData(3, fTrackCuts);
	196	PostData(4, fHelperPID);
	197	}
	198
	199	//________________________________________________________________________
	200
	201	void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
	202	{
	203	//Printf("An event");
	204	// main event loop
	205	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	206	if (!fAOD) {
	207	AliWarning("ERROR: AliAODEvent not available \n");
	208	return;
	209	}
86bac2e4	210
c683985a	211	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
86bac2e4	212	{
	213	AliFatal("Not processing AODs");
	214	}
	215
c683985a	216	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
	217
	218	//Default TPC priors
13dab2da	219	if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME we should modify the task to change priors
86bac2e4	220
ed33b9ef	221	Double_t Qvec=0.;
	222	if(fIsQvecCalibMode){
	223	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	224	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
a72cae45	225	else if (fVZEROside==2)Qvec=fEventCuts->GetqTPC();
ed33b9ef	226	}
ed33b9ef	227	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
86bac2e4	228
ed33b9ef	229	Double_t QvecMC = 0.;
ed33b9ef	230	if(fIsMC){
eba59d12	231	if(fIsQvecCalibMode){
725da720	232	QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside, fQgenType);
eba59d12	233	}
725da720	234	else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
c683985a	235	}
86bac2e4	236
509a25f9	237	Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
86bac2e4	238
c683985a	239	// First do MC to fill up the MC particle array
	240	TClonesArray *arrayMC = 0;
	241	if (fIsMC)
86bac2e4	242	{
	243	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	244	if (!arrayMC) {
	245	AliFatal("Error: MC particles branch not found!\n");
	246	}
	247	Int_t nMC = arrayMC->GetEntries();
	248	for (Int_t iMC = 0; iMC < nMC; iMC++)
c683985a	249	{
86bac2e4	250	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);
	251	if(!partMC->Charge()) continue;//Skip neutrals
	252	if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
	253
	254	//flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
	255	Double_t etaselected=-1.;
	256	if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
	257
	258	//pt cent IDgen isph y
	259	Double_t varSt[7];
	260	varSt[0]=partMC->Pt();
	261	varSt[1]=Cent;
	262	varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	263	varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	264	varSt[4]=partMC->Y();
	265	varSt[5]=etaselected;
	266
	267	if(fQst==0) varSt[6]=Qvec;
	268	else varSt[6]=QvecMC;
	269
	270	((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
	271
	272	//Printf("a particle");
	273
c683985a	274	}
86bac2e4	275	}
86bac2e4	276
c683985a	277	//main loop on tracks
86bac2e4	278
e19bc133	279	Int_t Nch = 0.;
86bac2e4	280
c683985a	281	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
f15c1f69	282	AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
f15c1f69	283	if(!track) AliFatal("Not a standard AOD");
c683985a	284	if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
c683985a	285	if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
dbbc07ac	286	if(!fFillOnlyEvents){
	287	Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
	288	Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
	289	Int_t IDgen=kSpUndefined;//set if MC
	290	Int_t isph=-999;
86bac2e4	291	// Int_t iswd=-999;
86bac2e4	292
dbbc07ac	293	if (arrayMC) {
86bac2e4	294	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));
	295	if (!partMC) {
	296	AliError("Cannot get MC particle");
	297	continue;
	298	}
	299	IDgen=fHelperPID->GetParticleSpecies(partMC);
	300	isph=partMC->IsPhysicalPrimary();
	301	//iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
	302
	303	if(fIsAOD160){// enabled for new ADO160 only
	304	if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
	305	if(partMC->IsSecondaryFromMaterial()) isph=3.;
	306
	307	}
c683985a	308	}
86bac2e4	309
93c24428	310	/* DCA */
93c24428	311	Double_t dcaxy = -999.;
3dd5ed48	312	Double_t dcaz = -999.;
86bac2e4	313
665f3db3	314	Double_t p[2];
3dd5ed48	315	if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
86bac2e4	316
3dd5ed48	317	if(dcaz >= fDCAzCut) continue;
896040ab	318
	319	//if the q vector is done using the TPC, we avoid overlap
	320	if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)continue;
86bac2e4	321
	322	//pt cent Q vec IDrec IDgen isph y
	323	// Double_t varTrk[8];
	324	Double_t varTrk[7];
dbbc07ac	325	varTrk[0]=track->Pt();
dbbc07ac	326	varTrk[1]=Cent;
09f5725e	327	if(fIsMC && fQtrk==1) varTrk[2]=QvecMC;
86bac2e4	328	else varTrk[2]=Qvec;
dbbc07ac	329	varTrk[3]=(Double_t)IDrec;
	330	varTrk[4]=(Double_t)IDgen;
	331	varTrk[5]=(Double_t)isph;
eba59d12	332	varTrk[6]=y;
86bac2e4	333	//varTrk[7]=dcaxy;
dbbc07ac	334	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
86bac2e4	335
dbbc07ac	336	//for nsigma PID fill double counting of ID
dbbc07ac	337	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
86bac2e4	338	Bool_t *HasDC;
	339	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	340	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	341	if(HasDC[ipart]==kTRUE){
	342	varTrk[3]=(Double_t)ipart;
	343	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	344	}
	345	}
c683985a	346	}
86bac2e4	347
dbbc07ac	348	//fill all charged (3)
	349	varTrk[3]=3.;
	350	varTrk[4]=3.;
	351	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	352	}//end if fFillOnlyEvents
86bac2e4	353
c683985a	354	//Printf("a track");
86bac2e4	355
dbbc07ac	356	Nch++;
c683985a	357	} // end loop on tracks
86bac2e4	358
ed33b9ef	359	Double_t varEv[4];
c683985a	360	varEv[0]=Cent;
c683985a	361	varEv[1]=Qvec;
e19bc133	362	varEv[2]=Nch;
c683985a	363	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
86bac2e4	364
c683985a	365	PostData(1, fOutput );
	366	PostData(2, fEventCuts);
	367	PostData(3, fTrackCuts);
	368	PostData(4, fHelperPID);
	369	}
	370
665f3db3	371	//_________________________________________________________________
665f3db3	372	Bool_t AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
86bac2e4	373
665f3db3	374	//AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
	375	//FIXME should update EventCuts?
	376	//FIXME add track->GetXYZ(p) method
86bac2e4	377
665f3db3	378	double xyz[3],cov[3];
86bac2e4	379
665f3db3	380	if (!trk->GetXYZ(xyz)) { // dca is not stored
	381	AliExternalTrackParam etp;
	382	etp.CopyFromVTrack(trk);
	383	AliVEvent* ev = (AliVEvent*)trk->GetEvent();
	384	if (!ev) {/printf("Event is not connected to the track\n");/ return kFALSE;}
	385	if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
	386	}
	387	p[0] = xyz[0];
	388	p[1] = xyz[1];
	389	return kTRUE;
	390
	391	}
	392
c683985a	393	//_________________________________________________________________
	394	void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
	395	{
	396	// Terminate
	397	}