Test task for Nano and standard AOD
authormfloris <michele.floris@cern.ch>
Tue, 18 Mar 2014 15:28:12 +0000 (16:28 +0100)
committerhristov <Peter.Hristov@cern.ch>
Thu, 27 Mar 2014 15:25:13 +0000 (16:25 +0100)
PWG/DevNanoAOD/AddTaskSpectraAllChNanoAOD.C [new file with mode: 0644]
PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.cxx [new file with mode: 0644]
PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.h [new file with mode: 0644]
PWG/DevNanoAOD/runAllCh.C [new file with mode: 0644]

diff --git a/PWG/DevNanoAOD/AddTaskSpectraAllChNanoAOD.C b/PWG/DevNanoAOD/AddTaskSpectraAllChNanoAOD.C
new file mode 100644 (file)
index 0000000..1d13ea6
--- /dev/null
@@ -0,0 +1,101 @@
+AliAnalysisTaskSpectraAllChNanoAOD* AddTaskSpectraAllChAOD(
+                                                      Bool_t mc=kFALSE,
+                                                      Double_t CentCutMin=0,
+                                                      Double_t CentCutMax=100,
+                                                      Double_t QvecCutMin=0,
+                                                      Double_t QvecCutMax=100,
+                                                      Double_t EtaMin=-0.8,
+                                                      Double_t EtaMax=0.8,
+                                                      Double_t pt=50.,
+                                                      Double_t ptTofMatch=.6,
+                                                       UInt_t trkbit=1,
+                                                      Double_t DCA=100000,
+                                                      UInt_t minNclsTPC=70,
+                                                      Double_t nsigmacut=5.,
+                                                      Int_t PIDtype=3,
+                                                      TString opt=""){
+  
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) 
+    {
+      ::Error("AddAliAnalysisTaskSpectraAllChAOD", "No analysis manager to connect to.");
+      return NULL;
+    }   
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) 
+    {
+      ::Error("AliAnalysisTaskSpectraAllChAOD", "This task requires an input event handler");
+      return NULL;
+    }   
+  
+  TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  if(type.Contains("ESD"))
+    {
+      ::Error("AliAnalysisTaskSpectraAllChAOD", "This task requires to run on AOD");
+      return NULL;
+    }
+  
+  AliSpectraAODTrackCuts  * trcuts = new AliSpectraAODTrackCuts(Form("TrackCuts%s",opt.Data()));  
+  trcuts->SetDCA(DCA);
+  trcuts->SetTrackBits(trkbit);
+  trcuts->SetPt(pt);
+  trcuts->SetPtTOFMatching(ptTofMatch);   
+  trcuts->SetEta(EtaMin,EtaMax);
+  trcuts->SetMinTPCcls(minNclsTPC);
+  trcuts->PrintCuts();
+  
+  AliSpectraAODEventCuts * evcuts = new AliSpectraAODEventCuts(Form("EventCuts%s",opt.Data()));
+  evcuts->SetQVectorCut(QvecCutMin,QvecCutMax);
+  evcuts->SetCentralityCutMax(CentCutMax);  
+  evcuts->SetCentralityCutMin(CentCutMin);
+  if(mc==1)evcuts->SetIsMC(kTRUE);
+  evcuts->PrintCuts();
+  
+  AliHelperPID *pid=new AliHelperPID();
+  pid->SetName(Form("HelperPID%s",opt.Data()));
+  pid->SetNSigmaCut(nsigmacut);
+  pid->SetPIDType(PIDtype);
+  if(PIDtype==3){
+    AliPIDCombined *pidc=new AliPIDCombined();
+    pidc->SetDefaultTPCPriors();
+    pid->SetPIDCombined(pidc);
+  }
+  
+  AliAnalysisTaskSpectraAllChNanoAOD *task = new AliAnalysisTaskSpectraAllChNanoAOD(Form("TaskAODSpectraCent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_TrBit%d%s",        
+                                                                                         CentCutMin,
+                                                                                         CentCutMax,
+                                                                                         QvecCutMin,
+                                                                                         QvecCutMax,
+                                                                                         EtaMin,
+                                                                                         EtaMax,
+                                                                                         trkbit,
+                                                                                         opt.Data()));
+  task->SetEventCuts(evcuts);
+  task->SetTrackCuts(trcuts);
+  task->SetHelperPID(pid);
+  if(mc==1)task->SetIsMC(kTRUE);
+  
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();
+  
+  TString typeofdata=mc?"MC":"Data";
+  
+  outputFileName += Form(":SpectraESE_%s%s",typeofdata.Data(),opt.Data());
+  
+  cout<<"outputFileName:  "<<outputFileName<<endl;
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();      
+  AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer(Form("chist%s",opt.Data()),                      TList::Class(),     AliAnalysisManager::kOutputContainer,outputFileName);
+  AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer(Form("cvcut%s",opt.Data()), AliSpectraAODEventCuts::Class(),     AliAnalysisManager::kOutputContainer,outputFileName);
+  AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer(Form("ctcut%s",opt.Data()), AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, outputFileName);
+  AliAnalysisDataContainer *coutputpt4 = mgr->CreateContainer(Form("cpid%s",opt.Data()),               AliHelperPID::Class(),     AliAnalysisManager::kOutputContainer, outputFileName);
+  mgr->AddTask(task);
+  
+  mgr->ConnectInput(task, 0, cinput);
+  mgr->ConnectOutput(task, 1, coutputpt1);
+  mgr->ConnectOutput(task, 2, coutputpt2);
+  mgr->ConnectOutput(task, 3, coutputpt3);
+  mgr->ConnectOutput(task, 4, coutputpt4);
+  
+  return task;
+}
diff --git a/PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.cxx b/PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.cxx
new file mode 100644 (file)
index 0000000..0a6019c
--- /dev/null
@@ -0,0 +1,371 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//         AliAnalysisTaskSpectraAllChNanoAOD class
+// NanoAOD edition: this is only meant to test the developement
+// version of NanoAODs
+//
+//-----------------------------------------------------------------
+
+#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 "AliAnalysisTaskSpectraAllChNanoAOD.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>
+#include "AliNanoAODHeader.h"
+#include "AliNanoAODTrack.h"
+
+using namespace AliHelperPIDNameSpace;
+using namespace std;
+
+ClassImp(AliAnalysisTaskSpectraAllChNanoAOD)
+
+//________________________________________________________________________
+AliAnalysisTaskSpectraAllChNanoAOD::AliAnalysisTaskSpectraAllChNanoAOD(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(40),
+  fnNchBins(200)
+{
+  // Default constructor
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, AliSpectraAODEventCuts::Class());
+  DefineOutput(3, AliSpectraAODTrackCuts::Class());
+  DefineOutput(4, AliHelperPID::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSpectraAllChNanoAOD::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=3;
+  //                                             cent             Q vec                Nch
+  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
+  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
+  Double_t xmaxHistRealEv[nvarev] = {       100.,               8.,               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 AliAnalysisTaskSpectraAllChNanoAOD::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");
+    } 
+
+  AliNanoAODHeader * headNano = dynamic_cast<AliNanoAODHeader*>((TObject*)fAOD->GetHeader());
+  
+  Bool_t isNano = (headNano != 0);
+  if(!isNano) {
+    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(!isNano) {
+      if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
+      else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
+    } else {
+      
+      if(fVZEROside==0)      Qvec=headNano->GetVar(0); // FIXME: we need proper getters here
+      else if (fVZEROside==1)Qvec=headNano->GetVar(1);
+      
+    }
+  }
+
+  Double_t Cent=isNano ? headNano->GetVar(2) : 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
+  
+  Int_t Nch = 0.;
+  
+  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
+    AliVTrack* track = (AliVTrack*) fAOD->GetTrack(iTracks);
+    if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge 
+    if(!isNano) {
+      if (!fTrackCuts->IsSelected((AliAODTrack*)track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
+    }
+    
+    if(!fFillOnlyEvents){
+      Int_t IDrec=isNano ? GetNanoTrackID (track) : fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector      
+      Double_t y= 0;
+      if(isNano) y = ((AliNanoAODTrack*)track)->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
+      else y = ((AliAODTrack*)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 (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[3];
+  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);
+}
+
+//_________________________________________________________________
+void   AliAnalysisTaskSpectraAllChNanoAOD::Terminate(Option_t *)
+{
+  // Terminate
+}
+
+
+Int_t AliAnalysisTaskSpectraAllChNanoAOD::GetNanoTrackID(AliVTrack * track) {
+  // Applies nsigma PID to nano tracks
+  AliNanoAODTrack * nanoTrack = dynamic_cast<AliNanoAODTrack*>(track);
+  if(!nanoTrack) AliFatal("Not a nano AOD track");
+
+  // Cache indexes
+  static const Int_t kcstNSigmaTPCPi  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPi");
+  static const Int_t kcstNSigmaTPCKa  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCKa");
+  static const Int_t kcstNSigmaTPCPr  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPr");
+  static const Int_t kcstNSigmaTOFPi  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPi");
+  static const Int_t kcstNSigmaTOFKa  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFKa");
+  static const Int_t kcstNSigmaTOFPr  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPr");
+
+  Double_t nSigmaPID = 3.0;
+
+
+
+  //get the identity of the particle with the minimum Nsigma
+  Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
+  if(nanoTrack->Pt() > fTrackCuts->GetPtTOFMatching()) {
+    nsigmaProton =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPr)*nanoTrack->GetVar(kcstNSigmaTPCPr)+nanoTrack->GetVar(kcstNSigmaTOFPr)*nanoTrack->GetVar(kcstNSigmaTOFPr));
+    nsigmaKaon   =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCKa)*nanoTrack->GetVar(kcstNSigmaTPCKa)+nanoTrack->GetVar(kcstNSigmaTOFKa)*nanoTrack->GetVar(kcstNSigmaTOFKa));
+    nsigmaPion   =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPi)*nanoTrack->GetVar(kcstNSigmaTPCPi)+nanoTrack->GetVar(kcstNSigmaTOFPi)*nanoTrack->GetVar(kcstNSigmaTOFPi));
+  }
+  else {
+    nsigmaProton =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPr));
+    nsigmaKaon   =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCKa));  
+    nsigmaPion   =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPi));  
+  }
+
+  
+
+  // guess the particle based on the smaller nsigma (within nSigmaPID)
+  if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
+  
+  if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < nSigmaPID)){
+    return kSpKaon;
+  }
+  if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < nSigmaPID)){    
+    return kSpPion;
+  }
+  if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < nSigmaPID)){
+    return kSpProton;
+  }
+  // else, return undefined
+  return kSpUndefined;
+  
+
+}
diff --git a/PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.h b/PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.h
new file mode 100644 (file)
index 0000000..fe5ca5d
--- /dev/null
@@ -0,0 +1,101 @@
+#ifndef ALIANALYSISTASKSPECTRAALLCHNANOAOD_H
+#define ALIANALYSISTASKSPECTRAALLCHNANOAOD_H
+
+/*  See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+//                      AliAnalysisTaskSpectraAllChNanoAOD
+//
+//
+//
+//
+// Author: Leonardo Milano, CERN
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class AliSpectraAODTrackCuts;
+class AliSpectraAODEventCuts;
+class AliHelperPID;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskSpectraAllChNanoAOD : public AliAnalysisTaskSE
+{
+ public:
+  // constructors
+ AliAnalysisTaskSpectraAllChNanoAOD() : AliAnalysisTaskSE(),
+    fAOD(0x0),
+    fTrackCuts(0x0),
+    fEventCuts(0x0),
+    fHelperPID(0x0),
+    fIsMC(0),
+    fDoDoubleCounting(0),
+    fFillOnlyEvents(0),
+    fCharge(0),
+    fVZEROside(0),
+    fOutput(0x0),
+    fnCentBins(20),
+    fnQvecBins(40),
+    fnNchBins(200)
+      {}
+  AliAnalysisTaskSpectraAllChNanoAOD(const char *name);
+  virtual ~AliAnalysisTaskSpectraAllChNanoAOD() {
+    Printf("calling detructor of AliAnalysisTaskSpectraAllChNanoAOD - To be implemented");
+  }
+  
+  void SetIsMC(Bool_t isMC = kFALSE)    {fIsMC = isMC; };
+  Bool_t GetIsMC()           const           { return fIsMC;};
+  void SetDoDoubleCounting(Bool_t doDoubleCounting = kFALSE)    {fDoDoubleCounting = doDoubleCounting; };
+  Bool_t GetDoDoubleCounting()           const           { return fDoDoubleCounting;};
+  void SetFillOnlyEvents(Bool_t fillOnlyEvents = kFALSE)    {fFillOnlyEvents = fillOnlyEvents; };
+  Bool_t GetFillOnlyEvents()           const           { return fFillOnlyEvents;};
+  void SetCharge(Int_t charge = 0)    {fCharge = charge; };
+  Int_t GetCharge()           const           { return fCharge;};
+  
+  void SetVZEROside(Int_t side = 0)    {fVZEROside = side; };
+  Int_t GetVZEROside()           const           { return fVZEROside;};
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  
+  AliSpectraAODTrackCuts      * GetTrackCuts()         {  return fTrackCuts; }
+  AliSpectraAODEventCuts      * GetEventCuts()         {  return fEventCuts; }
+  AliHelperPID                   * GetHelperPID()          { return fHelperPID; }
+  TList                          * GetOutputList()         { return fOutput; }
+  
+  void SetTrackCuts(AliSpectraAODTrackCuts * tc)       { fTrackCuts = tc; }
+  void SetEventCuts(AliSpectraAODEventCuts * vc)       { fEventCuts = vc; }
+  void SetHelperPID(AliHelperPID* pid)                     { fHelperPID = pid; }
+  void SetnCentBins(Int_t val)                             { fnCentBins = val; }
+  void SetnQvecBins(Int_t val)                             { fnQvecBins = val; }
+  void SetnNchBins(Int_t val)                             { fnNchBins = val; }
+
+
+  Int_t GetNanoTrackID(AliVTrack * track) ;
+  
+ private:
+  
+  AliAODEvent                   * fAOD;                         //! AOD object
+  AliSpectraAODTrackCuts      * fTrackCuts;                   // Track Cuts
+  AliSpectraAODEventCuts      * fEventCuts;                   // Event Cuts
+  AliHelperPID                   * fHelperPID;                    // points to class for PID
+  Bool_t                          fIsMC;                         // true if processing MC
+  Bool_t                          fDoDoubleCounting;           // true is double counting for Nsigma accepetd
+  Bool_t                          fFillOnlyEvents;               // if true fill only NSparseHistEv
+  Int_t                            fCharge;                      // charge to be selected
+  Int_t                            fVZEROside;                  // 0: VZERO-A 1: VZERO-C
+  TList                          * fOutput;                     // output list
+  Int_t                            fnCentBins;                  // number of bins for the centrality axis
+  Int_t                            fnQvecBins;                 // number of bins for the q vector axis
+  Int_t                            fnNchBins;                 // number of bins for the Nch axis
+  AliAnalysisTaskSpectraAllChNanoAOD(const AliAnalysisTaskSpectraAllChNanoAOD&);
+  AliAnalysisTaskSpectraAllChNanoAOD& operator=(const AliAnalysisTaskSpectraAllChNanoAOD&);
+  
+  ClassDef(AliAnalysisTaskSpectraAllChNanoAOD, 6);
+};
+
+#endif
diff --git a/PWG/DevNanoAOD/runAllCh.C b/PWG/DevNanoAOD/runAllCh.C
new file mode 100644 (file)
index 0000000..bb71279
--- /dev/null
@@ -0,0 +1,314 @@
+// run.C
+//
+// Template run macro for AliBasicTask.cxx/.h with example layout of
+// physics selections and options, in macro and task.
+//
+// Author: Arvinder Palaha
+//
+class AliAnalysisGrid;
+
+//______________________________________________________________________________
+void runAllCh(
+              Bool_t isNano = 1,
+              const char* runtype = "local", // local, proof or grid
+              const char *gridmode = "test", // Set the run mode (can be "full", "test", "offline", "submit" or "terminate"). Full & Test work for proof
+              const bool bMCtruth = 1, // 1 = MCEvent handler is on (MC truth), 0 = MCEvent handler is off (MC reconstructed/real data)
+              const bool bMCphyssel = 0, // 1 = looking at MC truth or reconstructed, 0 = looking at real data
+              const Long64_t nentries = 2000, // for local and proof mode, ignored in grid mode. Set to 1234567890 for all events.
+              const Long64_t firstentry = 0, // for local and proof mode, ignored in grid mode
+              const char *proofdataset = "/alice/data/LHC10c_000120821_p1", // path to dataset on proof cluster, for proof analysis
+              const char *proofcluster = "alice-caf.cern.ch", // which proof cluster to use in proof mode
+              const char *taskname = "example_task" // sets name of grid generated macros
+         )
+{
+    // check run type
+    if(runtype != "local" && runtype != "proof" && runtype != "grid"){
+        Printf("\n\tIncorrect run option, check first argument of run macro");
+        Printf("\tint runtype = local, proof or grid\n");
+        return;
+    }
+    Printf("%s analysis chosen",runtype);
+  
+    // load libraries
+  gSystem->Load("libCore.so");  
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libVMC");
+  gSystem->Load("libTree");
+  gSystem->Load("libProof");
+  gSystem->Load("libMatrix");
+  gSystem->Load("libMinuit");
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libESD");
+  gSystem->Load("libAOD");
+  gSystem->Load("libANALYSIS");
+  //  return;
+    gSystem->Load("libOADB");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libTENDER");
+  gSystem->Load("libCORRFW");
+  //  gSystem->Load("libPWG0base");
+  gSystem->Load("libMinuit");
+  gSystem->Load("libPWGTools");
+  gSystem->Load("libPWGLFspectra");
+  gSystem->Load("libPWGLFthermalfits");    
+  gSystem->Load("libPWGDevNanoAOD.so");
+
+  
+    // add aliroot indlude path
+    gROOT->ProcessLine(Form(".include %s/include",gSystem->ExpandPathName("$ALICE_ROOT")));
+    gROOT->SetStyle("Plain");
+        
+    // analysis manager
+    AliAnalysisManager* mgr = new AliAnalysisManager(taskname);
+    
+    // create the alien handler and attach it to the manager
+    AliAnalysisGrid *plugin = CreateAlienHandler(taskname, gridmode, proofcluster, proofdataset); 
+    mgr->SetGridHandler(plugin);
+    
+    AliInputEventHandler* iH = new AliAODInputHandler();
+    if(isNano) {
+      iH->SetEventSelection(new AliAnalysisNanoAODTrackCuts); // FIXME: we need this, otherwise we crash in AliAODInputHandler::BeginEvent where fIsSelectedResult = fEvent->GetHeader()->GetOfflineTrigger(). This is a temporary hack. In the future, it will be solved by using AliVHeader in the AliAODInputHandler.
+    }
+
+//    AliAODInputHandler* iH = new AliAODInputHandler();
+//    iH->SetInactiveBranches("tracks. vertices. v0s. cascades. jets. caloClusters. fmdClusters. pmdClusters. dimuons. AliAODZDC");
+//    iH->SetInactiveBranches("*");
+//    iH->SetCheckStatistics(kTRUE);
+    mgr->SetInputEventHandler(iH);
+        
+    // mc event handlerrunEx01.C
+    if(bMCtruth) {
+        AliMCEventHandler* mchandler = new AliMCEventHandler();
+        // Not reading track references
+        mchandler->SetReadTR(kFALSE);
+        mgr->SetMCtruthEventHandler(mchandler);
+    }   
+
+    // === Physics Selection Task ===
+    //
+    // In SelectCollisionCandidate(), default is kMB, so the task UserExec() 
+    // function is only called for these events.
+    // Options are:
+    //    kMB             Minimum Bias trigger
+    //    kMBNoTRD        Minimum bias trigger where the TRD is not read out
+    //    kMUON           Muon trigger
+    //    kHighMult       High-Multiplicity Trigger
+    //    kUserDefined    For manually defined trigger selection
+    //
+    // Multiple options possible with the standard AND/OR operators && and ||
+    // These all have the usual offline SPD or V0 selections performed.
+    //
+    // With a pointer to the physics selection object using physSelTask->GetPhysicsSelection(),
+    // one can manually set the selected and background classes using:
+    //    AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL")
+    //    AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
+    //
+    // One can also specify multiple classes at once, or require a class to NOT
+    // trigger, for e.g.
+    //    AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
+    //
+    // NOTE that manually setting the physics selection overrides the standard
+    // selection, so it must be done in completeness.
+    //
+    // ALTERNATIVELY, one can make the physics selection inside the task
+    // UserExec().
+    // For this case, comment out the task->SelectCol.... line, 
+    // and see AliBasicTask.cxx UserExec() function for details on this.
+
+//    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+//    AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(bMCphyssel);
+//    if(!physSelTask) { Printf("no physSelTask"); return; }
+    //AliPhysicsSelection *physSel = physSelTask->GetPhysicsSelection();
+    //physSel->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL");// #3119 #769");
+                
+    // PID RESPONSE
+    if(!isNano){
+      gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+      AliAnalysisTaskPIDResponse *taskPID=AddTaskPIDResponse(bMCtruth);
+      taskPID->SetUseTPCEtaCorrection(kTRUE); 
+      taskPID->SetUserDataRecoPass(2);
+    }
+
+
+    // create task
+    gROOT->LoadMacro("AddTaskSpectraAllChNanoAOD.C");
+    //load calibration object for LHC10h
+    AliAnalysisTaskSpectraAllChNanoAOD* task =     AddTaskSpectraAllChAOD(bMCtruth, // MC
+                                                                          0, // cent min
+                                                                          100, //cent max
+                                                                          -999999., // q min
+                                                                          999999., // q max
+                                                                          -0.8, // eta min
+                                                                          0.8, //  eta max
+                                                                          15.0, // pt max
+                                                                          0.6, // pt tof
+                                                                          1024, // trkbit
+                                                                          100000, // dca
+                                                                          70, // tpc cls
+                                                                          3.0, // nsigma
+                                                                          2, // pid type, 2 = nsigma TPC TOF
+                                                                          "_EvtOnly_F1024_VZEROA");
+    TFile * fCalib=TFile::Open("calibV0New.root");
+    task->GetEventCuts()->SetCalibFile(fCalib);
+    task->GetEventCuts()->SetIsLHC10h(kTRUE);
+    task->GetEventCuts()->SetCentralityMethod("V0M");
+    
+    task->SetVZEROside(0); //0: VZERO-A 1: VZERO-C
+    
+    task->SetFillOnlyEvents(0);
+    task->SetnCentBins(100);
+    task->SetnQvecBins(400);
+    task->SetnNchBins(1);
+
+    //    task->SelectCollisionCandidates(AliVEvent::kMB); // if physics selection performed in UserExec(), this line should be commented
+    mgr->AddTask(task);
+    
+    // set output root file name for different analysis
+    // TString outfilename = Form("list.%s.root",runtype);
+  
+    // create containers for input/output
+    // AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+    // AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("coutput1", TList::Class(), AliAnalysisManager::kOutputContainer, outfilename);
+        
+    // connect input/output
+    // mgr->ConnectInput(task, 0, cinput);
+    // mgr->ConnectOutput(task, 1, coutput1);
+        
+    // enable debug printouts
+    mgr->SetDebugLevel(10);
+    //    mgr->SetNSysInfo(100);
+    if (!mgr->InitAnalysis()) return;
+    mgr->PrintStatus();
+  
+    // start analysis
+    Printf("Starting Analysis....");
+    mgr->StartAnalysis(runtype,nentries,firstentry);
+}
+
+//______________________________________________________________________________
+AliAnalysisGrid* CreateAlienHandler(const char *taskname, const char *gridmode, const char *proofcluster, const char *proofdataset)
+{
+    AliAnalysisAlien *plugin = new AliAnalysisAlien();
+    // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+    plugin->SetRunMode(gridmode);
+
+    // Set versions of used packages
+    plugin->SetAPIVersion("V1.1x");
+    plugin->SetROOTVersion("v5-34-08-5");
+    plugin->SetAliROOTVersion("vAN-20140317");
+
+    // Declare input data to be processed.
+//    plugin->SetCheckCopy(kFALSE);
+
+    // Method 1: Create automatically XML collections using alien 'find' command.
+    // Define production directory LFN
+    plugin->SetGridDataDir("/alice/data/2010/LHC10b");
+    // On real reconstructed data:
+    // plugin->SetGridDataDir("/alice/data/2009/LHC09d");
+    // Set data search pattern
+    //plugin->SetDataPattern("*ESDs.root"); // THIS CHOOSES ALL PASSES
+    // Data pattern for reconstructed data
+    plugin->SetDataPattern("*ESDs/pass2/*ESDs.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH
+    //    plugin->SetDataPattern("ESDs/pass2/AOD038/*AliAOD.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH
+    plugin->SetRunPrefix("000");   // real data
+    // ...then add run numbers to be considered
+    Int_t runlist[15]={117039, 146859, 146858, 146856, 146824, 146817, 146806, 146805, 146804, 146803, 146802, 146801, 146748, 146747, 146746};  
+    for (Int_t ind=0; ind<1; ind++) {
+//     plugin->AddRunNumber(138275);
+      plugin->AddRunNumber(runlist[ind]);
+    }
+    //plugin->SetRunRange(114917,115322);
+    plugin->SetNrunsPerMaster(10); // 1
+    plugin->SetOutputToRunNo();
+    // comment out the next line when using the "terminate" option, unless
+    // you want separate merged files for each run
+    plugin->SetMergeViaJDL();
+
+    // Method 2: Declare existing data files (raw collections, xml collections, root file)
+    // If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir())
+    // XML collections added via this method can be combined with the first method if
+    // the content is compatible (using or not tags)
+    //   plugin->AddDataFile("tag.xml");
+    //   plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root");
+
+    // Define alien work directory where all files will be copied. Relative to alien $HOME.
+    plugin->SetGridWorkingDir(taskname);
+
+    // Declare alien output directory. Relative to working directory.
+    plugin->SetGridOutputDir("out"); // In this case will be $HOME/taskname/out
+
+    // Declare the analysis source files names separated by blancs. To be compiled runtime
+    // using ACLiC on the worker nodes.
+    plugin->SetAnalysisSource("AliAnalysisTaskEx01.cxx");
+
+    // Declare all libraries (other than the default ones for the framework. These will be
+    // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.
+    plugin->SetAdditionalLibs("AliAnalysisTaskEx01.h AliAnalysisTaskEx01.cxx");
+
+    // Declare the output file names separated by blancs.
+    // (can be like: file.root or file.root@ALICE::Niham::File)
+    // To only save certain files, use SetDefaultOutputs(kFALSE), and then
+    // SetOutputFiles("list.root other.filename") to choose which files to save
+    plugin->SetDefaultOutputs();
+    //plugin->SetOutputFiles("list.root");
+
+    // Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+    plugin->SetAnalysisMacro(Form("%s.C",taskname));
+
+    // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+    plugin->SetSplitMaxInputFileNumber(100);
+
+    // Optionally modify the executable name (default analysis.sh)
+    plugin->SetExecutable(Form("%s.sh",taskname));
+
+    // set number of test files to use in "test" mode
+    plugin->SetNtestFiles(10);
+
+    // Optionally resubmit threshold.
+    plugin->SetMasterResubmitThreshold(90);
+
+    // Optionally set time to live (default 30000 sec)
+    plugin->SetTTL(30000);
+
+    // Optionally set input format (default xml-single)
+    plugin->SetInputFormat("xml-single");
+
+    // Optionally modify the name of the generated JDL (default analysis.jdl)
+    plugin->SetJDLName(Form("%s.jdl",taskname));
+
+    // Optionally modify job price (default 1)
+    plugin->SetPrice(1);      
+
+    // Optionally modify split mode (default 'se')    
+    plugin->SetSplitMode("se");
+    
+    //----------------------------------------------------------
+    //---      PROOF MODE SPECIFIC SETTINGS         ------------
+    //---------------------------------------------------------- 
+    // Proof cluster
+    plugin->SetProofCluster(proofcluster);
+    // Dataset to be used   
+    plugin->SetProofDataSet(proofdataset);
+    // May need to reset proof. Supported modes: 0-no reset, 1-soft, 2-hard
+    plugin->SetProofReset(0);
+    // May limit number of workers
+    plugin->SetNproofWorkers(0);
+    // May limit the number of workers per slave
+    plugin->SetNproofWorkersPerSlave(1);   
+    // May use a specific version of root installed in proof
+    plugin->SetRootVersionForProof("current");
+    // May set the aliroot mode. Check http://aaf.cern.ch/node/83 
+    plugin->SetAliRootMode("default"); // Loads AF libs by default
+    // May request ClearPackages (individual ClearPackage not supported)
+    plugin->SetClearPackages(kFALSE);
+    // Plugin test mode works only providing a file containing test file locations, used in "local" mode also
+    plugin->SetFileForTestMode("files.txt"); // file should contain path name to a local directory containg *ESDs.root etc
+    // Request connection to alien upon connection to grid
+    plugin->SetProofConnectGrid(kFALSE);
+    // Other PROOF specific parameters
+    plugin->SetProofParameter("PROOF_UseMergers","-1");
+    printf("Using: PROOF_UseMergers   : %s\n", plugin->GetProofParameter("PROOF_UseMergers"));
+    return plugin;
+}
+