]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx
flag for double counting
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx
index 2ff52b05aa1e6c73d1f3def701ee11b68fb03e4d..3b89a83c5062b69417aab859c32b5f10760f0fac 100644 (file)
-/**************************************************************************\r
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-\r
-//-----------------------------------------------------------------\r
-//         AliAnalysisTaskSpectraAllChAOD class\r
-//-----------------------------------------------------------------\r
-\r
-#include "TChain.h"\r
-#include "TTree.h"\r
-#include "TLegend.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "THnSparse.h"\r
-#include "TCanvas.h"\r
-#include "AliAnalysisTask.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliVParticle.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODInputHandler.h"\r
-#include "AliAnalysisTaskSpectraAllChAOD.h"\r
-#include "AliAnalysisTaskESDfilter.h"\r
-#include "AliAnalysisDataContainer.h"\r
-#include "AliSpectraAODTrackCuts.h"\r
-#include "AliSpectraAODEventCuts.h"\r
-#include "AliHelperPID.h"\r
-#include "AliCentrality.h"\r
-#include "TProof.h"\r
-#include "AliVEvent.h"\r
-#include "AliStack.h"\r
-#include <TMCProcess.h>\r
-\r
-#include <iostream>\r
-\r
-using namespace AliHelperPIDNameSpace;\r
-using namespace std;\r
-\r
-ClassImp(AliAnalysisTaskSpectraAllChAOD)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  \r
-  fAOD(0x0),\r
-  fTrackCuts(0x0),\r
-  fEventCuts(0x0),\r
-  fHelperPID(0x0),\r
-  fIsMC(0),\r
-  fOutput(0x0)\r
-{\r
-  // Default constructor\r
-  DefineInput(0, TChain::Class());\r
-  DefineOutput(1, TList::Class());\r
-  DefineOutput(2, AliSpectraAODEventCuts::Class());\r
-  DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
-  DefineOutput(4, AliHelperPID::Class());\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()\r
-{\r
-  // create output objects\r
-  fOutput=new TList();\r
-  fOutput->SetOwner();\r
-  fOutput->SetName("fOutput");\r
-  \r
-  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
-  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
-  if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");\r
-  \r
-  // binning common to all the THn\r
-  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.};\r
-  const Int_t nptBins=26;\r
-  const Int_t nCentBins=20;\r
-  const Int_t nQvecBins=50;\r
-  \r
-  //dimensions of THnSparse for tracks\r
-  const Int_t nvartrk=8;\r
-  //                                             pt          cent        Q vec     IDrec     IDgen       isph           iswd      y\r
-  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, nCentBins,   nQvecBins,        3,        3,         2,          2,       2};\r
-  Double_t xminHistRealTrk[nvartrk] = {         0.,         0.,            0.,     -0.5,      -0.5,      -0.5,        -0.5,   -0.5};\r
-  Double_t xmaxHistRealTrk[nvartrk] = {       10.,      100.,          10.,      2.5,      2.5,       1.5,         1.5,     0.5};    \r
-  THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);\r
-  NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");\r
-  NSparseHistTrk->SetBinEdges(0,ptBins);\r
-  NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
-  NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");\r
-  NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");\r
-  NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");\r
-  NSparseHistTrk->GetAxis(5)->SetTitle("isph");\r
-  NSparseHistTrk->GetAxis(6)->SetTitle("iswd");\r
-  NSparseHistTrk->GetAxis(7)->SetTitle("y");\r
-  fOutput->Add(NSparseHistTrk);\r
-  \r
-  //dimensions of THnSparse for stack\r
-  const Int_t nvarst=5;\r
-  //                                             pt          cent    IDgen        isph        y\r
-  Int_t    binsHistRealSt[nvarst] = {      nptBins,   nCentBins,        3,         2,        2};\r
-  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5};\r
-  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5};\r
-  THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);\r
-  NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");\r
-  NSparseHistSt->SetBinEdges(0,ptBins);\r
-  NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
-  NSparseHistSt->GetAxis(2)->SetTitle("ID gen");\r
-  NSparseHistSt->GetAxis(3)->SetTitle("isph");\r
-  NSparseHistSt->GetAxis(4)->SetTitle("y");\r
-  fOutput->Add(NSparseHistSt);\r
-  \r
-  //dimensions of THnSparse for the normalization\r
-  const Int_t nvarev=2;\r
-  //                                             cent             Q vec   \r
-  Int_t    binsHistRealEv[nvarev] = {     nCentBins,       nQvecBins};\r
-  Double_t xminHistRealEv[nvarev] = {           0.,               0.};\r
-  Double_t xmaxHistRealEv[nvarev] = {       100.,              10.};\r
-  THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);\r
-  NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
-  NSparseHistEv->GetAxis(1)->SetTitle("Q vec");\r
-  fOutput->Add(NSparseHistEv);\r
-  \r
-  PostData(1, fOutput  );\r
-  PostData(2, fEventCuts);\r
-  PostData(3, fTrackCuts);\r
-  PostData(4, fHelperPID);\r
-}\r
-\r
-//________________________________________________________________________\r
-\r
-void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)\r
-{\r
-  //Printf("An event");\r
-  // main event loop\r
-  fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
-  if (!fAOD) {\r
-    AliWarning("ERROR: AliAODEvent not available \n");\r
-    return;\r
-  }\r
-  \r
-  if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
-    {\r
-      AliFatal("Not processing AODs");\r
-    }\r
-  \r
-  if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
-  \r
-  Double_t Qvec=0.;//in case of MC we save space in the memory\r
-  if(!fIsMC)Qvec=fEventCuts->GetqV0A();//FIXME we'll have to combine A and C\r
-  Double_t Cent=fEventCuts->GetCent();\r
-    \r
-  // First do MC to fill up the MC particle array\r
-  TClonesArray *arrayMC = 0;\r
-  if (fIsMC)\r
-    {\r
-      arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
-      if (!arrayMC) {\r
-       AliFatal("Error: MC particles branch not found!\n");\r
-      }\r
-      Int_t nMC = arrayMC->GetEntries();\r
-      for (Int_t iMC = 0; iMC < nMC; iMC++)\r
-       {\r
-         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
-         if(!partMC->Charge()) continue;//Skip neutrals\r
-         if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!\r
-         \r
-         //pt     cent    IDgen        isph        y\r
-         Double_t varSt[5];\r
-         varSt[0]=partMC->Pt();\r
-         varSt[1]=Cent;\r
-         varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);\r
-         varSt[3]=(Double_t)partMC->IsPhysicalPrimary();\r
-         varSt[4]=partMC->Y();\r
-         ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop\r
-         \r
-         //Printf("a particle");\r
-         \r
-       }\r
-    }\r
-  \r
-  //main loop on tracks\r
-  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {\r
-    AliAODTrack* track = fAOD->GetTrack(iTracks);\r
-    if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)\r
-    Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector\r
-    Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));\r
-    Int_t IDgen=kSpUndefined;//set if MC\r
-    Int_t isph=-999;\r
-    Int_t iswd=-999;\r
-    \r
-    if (arrayMC) {\r
-      AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
-      if (!partMC) { \r
-       AliError("Cannot get MC particle");\r
-       continue; \r
-      }\r
-      IDgen=fHelperPID->GetParticleSpecies(partMC);\r
-      isph=partMC->IsPhysicalPrimary();\r
-      iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions\r
-    }\r
-    \r
-    //pt     cent    Q vec     IDrec     IDgen       isph           iswd      y\r
-    Double_t varTrk[8];\r
-    varTrk[0]=track->Pt();\r
-    varTrk[1]=Cent;\r
-    varTrk[2]=Qvec;\r
-    varTrk[3]=(Double_t)IDrec;\r
-    varTrk[4]=(Double_t)IDgen;\r
-    varTrk[5]=(Double_t)isph;\r
-    varTrk[6]=(Double_t)iswd;\r
-    varTrk[7]=y;\r
-    ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
-    \r
-    //Printf("a track");\r
-    \r
-  } // end loop on tracks\r
-  \r
-  Double_t varEv[2];\r
-  varEv[0]=Cent;\r
-  varEv[1]=Qvec;\r
-  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop\r
-  \r
-  PostData(1, fOutput  );\r
-  PostData(2, fEventCuts);\r
-  PostData(3, fTrackCuts);\r
-  PostData(4, fHelperPID);\r
-}\r
-\r
-//_________________________________________________________________\r
-void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)\r
-{\r
-  // Terminate\r
-}\r
+/**************************************************************************
+ * 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),
+  fCharge(0),
+  fVZEROside(0),
+  fOutput(0x0),
+  fnCentBins(20),
+  fnQvecBins(40)
+{
+  // 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.};
+  const Int_t nptBins=26;
+  
+  //dimensions of THnSparse for tracks
+  const Int_t nvartrk=8;
+  //                                             pt          cent        Q vec     IDrec     IDgen       isph           iswd      y
+  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,   fnQvecBins,        4,        3,         2,          2,       2};
+  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,            0.,      -.5,      -0.5,      -0.5,        -0.5,   -0.5};
+  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,            8.,      3.5,      2.5,       1.5,         1.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("iswd");
+  NSparseHistTrk->GetAxis(6)->SetName("iswd");
+  NSparseHistTrk->GetAxis(7)->SetTitle("y");
+  NSparseHistTrk->GetAxis(7)->SetName("y");
+  fOutput->Add(NSparseHistTrk);
+  
+  //dimensions of THnSparse for stack
+  const Int_t nvarst=5;
+  //                                             pt          cent    IDgen        isph        y
+  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2};
+  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5};
+  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5};
+  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");
+  fOutput->Add(NSparseHistSt);
+  
+  //dimensions of THnSparse for the normalization
+  const Int_t nvarev=2;
+  //                                             cent             Q vec   
+  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins};
+  Double_t xminHistRealEv[nvarev] = {           0.,               0.};
+  Double_t xmaxHistRealEv[nvarev] = {       100.,              10.};
+  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");
+  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 maybe this can go in the UserCreateOutputObject?
+  
+  Double_t Qvec=0.;//in case of MC we save space in the memory
+  if(!fIsMC){
+    if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
+    else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
+  }
+
+  Double_t Cent=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
+         
+         if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
+         
+         //pt     cent    IDgen        isph        y
+         Double_t varSt[5];
+         varSt[0]=partMC->Pt();
+         varSt[1]=Cent;
+         varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
+         varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
+         varSt[4]=partMC->Y();
+         ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
+         
+         //Printf("a particle");
+         
+       }
+    }
+  
+  //main loop on tracks
+  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
+    AliAODTrack* track = fAOD->GetTrack(iTracks);
+    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)
+    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
+    }
+    
+    //pt     cent    Q vec     IDrec     IDgen       isph           iswd      y
+    Double_t varTrk[8];
+    varTrk[0]=track->Pt();
+    varTrk[1]=Cent;
+    varTrk[2]=Qvec;
+    varTrk[3]=(Double_t)IDrec;
+    varTrk[4]=(Double_t)IDgen;
+    varTrk[5]=(Double_t)isph;
+    varTrk[6]=(Double_t)iswd;
+    varTrk[7]=y;
+    ((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 (4)
+    varTrk[3]=4.;
+    ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
+    
+    //Printf("a track");
+    
+  } // end loop on tracks
+  
+  Double_t varEv[2];
+  varEv[0]=Cent;
+  varEv[1]=Qvec;
+  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
+  
+  PostData(1, fOutput  );
+  PostData(2, fEventCuts);
+  PostData(3, fTrackCuts);
+  PostData(4, fHelperPID);
+}
+
+//_________________________________________________________________
+void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
+{
+  // Terminate
+}