[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraAODEventCuts.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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliSpectraAODEventCuts class
//-----------------------------------------------------------------

#include "TChain.h"
#include "TTree.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TH1I.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisDataContainer.h"
#include "AliSpectraAODEventCuts.h"
#include "AliSpectraAODTrackCuts.h"
//#include "AliSpectraAODHistoManager.h"
#include <iostream>

using namespace std;

ClassImp(AliSpectraAODEventCuts)

AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0), fTrackBits(0),fIsMC(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
,fHistoEP(0)
{
  // Constructor
  fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
  fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",500,-15,15);
  fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",500,-15,15);
  fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
  fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
  fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
  //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
  //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
  fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
  fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
  fCentralityCutMin = 0.0;      // default value of centrality cut minimum, 0 ~ no cut
  fCentralityCutMax = 10000.0;  // default value of centrality cut maximum,  ~ no cut
  // fQVectorPosCutMin=0.0;
  // fQVectorPosCutMax=10000.0;
  // fQVectorNegCutMin=0.0;
  // fQVectorNegCutMax=10000.0;
  fQVectorCutMin=0.0;
  fQVectorCutMax=10000.0;
  fTrackBits=1;
}

//______________________________________________________
Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
{
  // Returns true if Event Cuts are selected and applied
  fAOD = aod;
  fTrackCuts = trackcuts;
  fHistoCuts->Fill(kProcessedEvents);
  Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
  if(!IsPhysSel)return IsPhysSel;
  fHistoCuts->Fill(kPhysSelEvents);
  //loop on tracks, before event selection, filling QA histos
  AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
  if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
  fIsSelected =kFALSE;
  if(CheckVtxRange() && CheckCentralityCut()){ //selection on vertex and Centrality
    fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
  }
  if(fIsSelected){
    fHistoCuts->Fill(kAcceptedEvents);
    if(vertex)fHistoVtxAftSel->Fill(vertex->GetZ());
  }
  Int_t Nch=0;
  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
    AliAODTrack* track = fAOD->GetTrack(iTracks);
    if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
    fHistoEtaBefSel->Fill(track->Eta());
    if(fIsSelected){
      fHistoEtaAftSel->Fill(track->Eta());
      Nch++;
    }
  }
  if(fIsSelected)fHistoNChAftSel->Fill(Nch);
  return fIsSelected;
}

//______________________________________________________
Bool_t AliSpectraAODEventCuts::CheckVtxRange()
{
  // reject events outside of range
  AliAODVertex * vertex = fAOD->GetPrimaryVertex();
  if (!vertex)
    {
      fHistoCuts->Fill(kVtxNoEvent);
      return kFALSE;
    }
  if (TMath::Abs(vertex->GetZ()) < 10)
    {
      return kTRUE;
    }
  fHistoCuts->Fill(kVtxRange);
  return kFALSE;
}

//______________________________________________________
Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
{
  // Check centrality cut
  Double_t cent=0;
  if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile("V0M");
  else cent=ApplyCentralityPatchAOD049();
  
  if ( (cent <= fCentralityCutMax)  &&  (cent >= fCentralityCutMin) )  return kTRUE;   
  fHistoCuts->Fill(kVtxCentral);
  return kFALSE;
}

//______________________________________________________
Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
{
  // Check qvector cut
  /// FIXME: Q vector
  // //Selection on QVector, before ANY other selection on the event
  // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
  // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
  // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
 
  // Int_t multPos = 0;
  // Int_t multNeg = 0;
  // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
  //   AliAODTrack* aodTrack = fAOD->GetTrack(iT);
  //   if (!aodTrack->TestFilterBit(fTrackBits)) continue;
  //   if (aodTrack->Eta() >= 0){
  //     multPos++;
  //     Qx2EtaPos += TMath::Cos(2*aodTrack->Phi()); 
  //     Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
  //   } else {
  //     multNeg++;
  //     Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi()); 
  //     Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
  //   }
  // } 
  // Double_t qPos=-999;
  // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
  // Double_t qNeg=-999;
  // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
  //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
  
  Double_t qxEPVZERO = 0, qyEPVZERO = 0;
  Double_t qVZERO = -999;
  Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
  
  qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
  if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
  fHistoQVector->Fill(qVZERO);
  fHistoEP->Fill(psi);
  
  fHistoCuts->Fill(kQVector);
  // fHistoQVectorPos->Fill(qPos);
  // fHistoQVectorNeg->Fill(qNeg);
  return kTRUE;
}

//______________________________________________________
void AliSpectraAODEventCuts::PrintCuts()
{
  // print info about event cuts
  cout << "Event Stats" << endl;
  cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
  cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
  cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
  cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
  cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
  cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
  cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
  cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
  // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
  // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
  cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
  cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
}
//______________________________________________________

Double_t AliSpectraAODEventCuts::ApplyCentralityPatchAOD049()
{
   //
   //Apply centrality patch for AOD049 outliers
   //
  // if (fCentralityType!="V0M") {
  //   AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
  //   return -999.0;
  // }
  AliCentrality *centrality = fAOD->GetCentrality();
   if (!centrality) {
     AliWarning("Cannot get centrality from AOD event.");
     return -999.0;
   }
   
   Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
   /*
     Bool_t isSelRun = kFALSE;
     Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
     if(cent<0){
     Int_t quality = centrality->GetQuality();
     if(quality<=1){
     cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
     } else {
     Int_t runnum=aodEvent->GetRunNumber();
     for(Int_t ir=0;ir<5;ir++){
     if(runnum==selRun[ir]){
     isSelRun=kTRUE;
     break;
     }
     }
     if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
     }
     }
   */
   if(cent>=0.0) {
     Float_t v0 = 0.0;
     AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
     AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
     v0+=aodV0->GetMTotV0A();
     v0+=aodV0->GetMTotV0C();
     if ( (cent==0) && (v0<19500) ) {
       AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
       return -999.0;
     }
     Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
     Float_t val = 1.30552 +  0.147931 * v0;
     
     Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
			      120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
			      92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
			      68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
			      51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
			      37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
			      26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
			      19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
			      12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
			      13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
     };
     
     if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
       AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
       return -999.0;
     }
   } else {
     //force it to be -999. whatever the negative value was
     cent = -999.;
   }
   return cent;
}

//______________________________________________________


Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
{
  // Merge a list of AliSpectraAODEventCuts objects with this.
  // Returns the number of merged objects (including this).

  //  AliInfo("Merging");

  if (!list)
    return 0;

  if (list->IsEmpty())
    return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj;

  // collections of all histograms
  TList collections;//FIXME we should use only 1 collection
  TList collections_histoVtxBefSel;
  TList collections_histoVtxAftSel;
  TList collections_histoEtaBefSel;
  TList collections_histoEtaAftSel;
  TList collections_histoNChAftSel;
  // TList collections_histoQVectorPos;
  // TList collections_histoQVectorNeg;
  TList collections_histoQVector;
  TList collections_histoEP;

  Int_t count = 0;

  while ((obj = iter->Next())) {
    AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
    if (entry == 0) 
      continue;

    TH1I * histo = entry->GetHistoCuts();      
    collections.Add(histo);
    TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();      
    collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
    TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();      
    collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
    TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();      
    collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
    TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();      
    collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
    TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();      
    collections_histoNChAftSel.Add(histo_histoNChAftSel);
    // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();      
    // collections_histoQVectorPos.Add(histo_histoQVectorPos);
    // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();      
    // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
    TH1F * histo_histoQVector = entry->GetHistoQVector();      
    collections_histoQVector.Add(histo_histoQVector);
    TH1F * histo_histoEP = entry->GetHistoEP();      
    collections_histoEP.Add(histo_histoEP);
    count++;
  }
  
  fHistoCuts->Merge(&collections);
  fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
  fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
  fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
  fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
  fHistoNChAftSel->Merge(&collections_histoNChAftSel);
  // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
  // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
  fHistoQVector->Merge(&collections_histoQVector);
  fHistoEP->Merge(&collections_histoEP);
  

  delete iter;

  return count+1;
}
Commit	Line	Data
c88234ad	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	//-----------------------------------------------------------------
	18	// AliSpectraAODEventCuts class
	19	//-----------------------------------------------------------------
	20
	21	#include "TChain.h"
	22	#include "TTree.h"
	23	#include "TLegend.h"
	24	#include "TH1F.h"
829b5a81	25	#include "TH1I.h"
c88234ad	26	#include "TH2F.h"
	27	#include "TCanvas.h"
	28	#include "AliAnalysisTask.h"
	29	#include "AliAnalysisManager.h"
	30	#include "AliAODTrack.h"
	31	#include "AliAODMCParticle.h"
	32	#include "AliAODEvent.h"
	33	#include "AliAODInputHandler.h"
	34	#include "AliAnalysisTaskESDfilter.h"
	35	#include "AliAnalysisDataContainer.h"
	36	#include "AliSpectraAODEventCuts.h"
f6a38178	37	#include "AliSpectraAODTrackCuts.h"
829b5a81	38	//#include "AliSpectraAODHistoManager.h"
c88234ad	39	#include <iostream>
	40
	41	using namespace std;
	42
	43	ClassImp(AliSpectraAODEventCuts)
	44
16ee3540	45	AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0), fTrackBits(0),fIsMC(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
16ee3540	46	,fHistoEP(0)
c88234ad	47	{
f6a38178	48	// Constructor
	49	fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
	50	fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",500,-15,15);
	51	fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",500,-15,15);
	52	fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
	53	fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
00493191	54	fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
16ee3540	55	//fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
	56	//fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
	57	fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
	58	fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
f6a38178	59	fCentralityCutMin = 0.0; // default value of centrality cut minimum, 0 ~ no cut
f6a38178	60	fCentralityCutMax = 10000.0; // default value of centrality cut maximum, ~ no cut
16ee3540	61	// fQVectorPosCutMin=0.0;
	62	// fQVectorPosCutMax=10000.0;
	63	// fQVectorNegCutMin=0.0;
	64	// fQVectorNegCutMax=10000.0;
	65	fQVectorCutMin=0.0;
	66	fQVectorCutMax=10000.0;
829b5a81	67	fTrackBits=1;
c88234ad	68	}
	69
	70	//______________________________________________________
f6a38178	71	Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts *trackcuts)
c88234ad	72	{
f6a38178	73	// Returns true if Event Cuts are selected and applied
	74	fAOD = aod;
	75	fTrackCuts = trackcuts;
	76	fHistoCuts->Fill(kProcessedEvents);
decf69d9	77	Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
decf69d9	78	if(!IsPhysSel)return IsPhysSel;
8961bb9a	79	fHistoCuts->Fill(kPhysSelEvents);
f6a38178	80	//loop on tracks, before event selection, filling QA histos
	81	AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
	82	if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
decf69d9	83	fIsSelected =kFALSE;
decf69d9	84	if(CheckVtxRange() && CheckCentralityCut()){ //selection on vertex and Centrality
00493191	85	fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
decf69d9	86	}
f6a38178	87	if(fIsSelected){
	88	fHistoCuts->Fill(kAcceptedEvents);
	89	if(vertex)fHistoVtxAftSel->Fill(vertex->GetZ());
	90	}
ae0fdd7d	91	Int_t Nch=0;
f6a38178	92	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
f6a38178	93	AliAODTrack* track = fAOD->GetTrack(iTracks);
decf69d9	94	if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
f6a38178	95	fHistoEtaBefSel->Fill(track->Eta());
ae0fdd7d	96	if(fIsSelected){
	97	fHistoEtaAftSel->Fill(track->Eta());
	98	Nch++;
	99	}
f6a38178	100	}
ae0fdd7d	101	if(fIsSelected)fHistoNChAftSel->Fill(Nch);
f6a38178	102	return fIsSelected;
c88234ad	103	}
	104
	105	//______________________________________________________
	106	Bool_t AliSpectraAODEventCuts::CheckVtxRange()
	107	{
	108	// reject events outside of range
f6a38178	109	AliAODVertex * vertex = fAOD->GetPrimaryVertex();
	110	if (!vertex)
	111	{
c88234ad	112	fHistoCuts->Fill(kVtxNoEvent);
c88234ad	113	return kFALSE;
f6a38178	114	}
	115	if (TMath::Abs(vertex->GetZ()) < 10)
	116	{
c88234ad	117	return kTRUE;
f6a38178	118	}
	119	fHistoCuts->Fill(kVtxRange);
	120	return kFALSE;
c88234ad	121	}
	122
	123	//______________________________________________________
	124	Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
	125	{
f6a38178	126	// Check centrality cut
972a21ad	127	Double_t cent=0;
936e2842	128	if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile("V0M");
972a21ad	129	else cent=ApplyCentralityPatchAOD049();
	130
	131	if ( (cent <= fCentralityCutMax) && (cent >= fCentralityCutMin) ) return kTRUE;
f6a38178	132	fHistoCuts->Fill(kVtxCentral);
f6a38178	133	return kFALSE;
c88234ad	134	}
c88234ad	135
decf69d9	136	//______________________________________________________
	137	Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
	138	{
	139	// Check qvector cut
	140	/// FIXME: Q vector
	141	// //Selection on QVector, before ANY other selection on the event
	142	// //Spectra MUST be normalized wrt events AFTER the selection on Qvector
16ee3540	143	// Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
	144	// Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
	145
	146	// Int_t multPos = 0;
	147	// Int_t multNeg = 0;
	148	// for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
	149	// AliAODTrack* aodTrack = fAOD->GetTrack(iT);
	150	// if (!aodTrack->TestFilterBit(fTrackBits)) continue;
	151	// if (aodTrack->Eta() >= 0){
	152	// multPos++;
	153	// Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
	154	// Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
	155	// } else {
	156	// multNeg++;
	157	// Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
	158	// Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
	159	// }
	160	// }
	161	// Double_t qPos=-999;
	162	// if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPosQx2EtaPos + Qy2EtaPosQy2EtaPos)/multPos);
	163	// Double_t qNeg=-999;
	164	// if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNegQx2EtaNeg + Qy2EtaNegQy2EtaNeg)/multNeg);
	165	//if(qPos<fQVectorPosCutMin \|\| qPos>fQVectorPosCutMax \|\| qNeg<fQVectorNegCutMin \|\| qNeg>fQVectorNegCutMax)return kFALSE;
decf69d9	166
16ee3540	167	Double_t qxEPVZERO = 0, qyEPVZERO = 0;
	168	Double_t qVZERO = -999;
	169	Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
00493191	170
16ee3540	171	qVZERO= TMath::Sqrt(qxEPVZEROqxEPVZERO + qyEPVZEROqyEPVZERO);
	172	if(qVZERO<fQVectorCutMin \|\| qVZERO>fQVectorCutMax)return kFALSE;
	173	fHistoQVector->Fill(qVZERO);
	174	fHistoEP->Fill(psi);
00493191	175
decf69d9	176	fHistoCuts->Fill(kQVector);
16ee3540	177	// fHistoQVectorPos->Fill(qPos);
16ee3540	178	// fHistoQVectorNeg->Fill(qNeg);
decf69d9	179	return kTRUE;
decf69d9	180	}
decf69d9	181
c88234ad	182	//______________________________________________________
	183	void AliSpectraAODEventCuts::PrintCuts()
	184	{
f6a38178	185	// print info about event cuts
	186	cout << "Event Stats" << endl;
	187	cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
	188	cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
936e2842	189	cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
f6a38178	190	cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
	191	cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
	192	cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
decf69d9	193	cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
829b5a81	194	cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
16ee3540	195	// cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
	196	// cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
	197	cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
829b5a81	198	cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
f6a38178	199	}
c88234ad	200	//______________________________________________________
c88234ad	201
972a21ad	202	Double_t AliSpectraAODEventCuts::ApplyCentralityPatchAOD049()
	203	{
	204	//
	205	//Apply centrality patch for AOD049 outliers
	206	//
	207	// if (fCentralityType!="V0M") {
	208	// AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
	209	// return -999.0;
	210	// }
	211	AliCentrality *centrality = fAOD->GetCentrality();
	212	if (!centrality) {
	213	AliWarning("Cannot get centrality from AOD event.");
	214	return -999.0;
	215	}
	216
	217	Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
	218	/*
	219	Bool_t isSelRun = kFALSE;
	220	Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
	221	if(cent<0){
	222	Int_t quality = centrality->GetQuality();
	223	if(quality<=1){
	224	cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
	225	} else {
	226	Int_t runnum=aodEvent->GetRunNumber();
	227	for(Int_t ir=0;ir<5;ir++){
	228	if(runnum==selRun[ir]){
	229	isSelRun=kTRUE;
	230	break;
	231	}
	232	}
	233	if((quality==8\|\|quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
	234	}
	235	}
	236	*/
	237	if(cent>=0.0) {
	238	Float_t v0 = 0.0;
	239	AliAODEvent aodEvent = (AliAODEvent )fAOD;
	240	AliAODVZERO aodV0 = (AliAODVZERO ) aodEvent->GetVZEROData();
	241	v0+=aodV0->GetMTotV0A();
	242	v0+=aodV0->GetMTotV0C();
	243	if ( (cent==0) && (v0<19500) ) {
	244	AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
	245	return -999.0;
	246	}
	247	Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
	248	Float_t val = 1.30552 + 0.147931 * v0;
	249
	250	Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
	251	120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
	252	92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
	253	68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
	254	51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
	255	37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
	256	26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
	257	19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
	258	12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
	259	13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
	260	};
	261
	262	if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
	263	AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
	264	return -999.0;
	265	}
266	} else {
267	//force it to be -999. whatever the negative value was
268	cent = -999.;
269	}
270	return cent;
271	}
272
273	//______________________________________________________
274
275
c88234ad	276	Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
	277	{
	278	// Merge a list of AliSpectraAODEventCuts objects with this.
	279	// Returns the number of merged objects (including this).
	280
	281	// AliInfo("Merging");
	282
	283	if (!list)
	284	return 0;
	285
	286	if (list->IsEmpty())
	287	return 1;
	288
	289	TIterator* iter = list->MakeIterator();
	290	TObject* obj;
	291
	292	// collections of all histograms
ae0fdd7d	293	TList collections;//FIXME we should use only 1 collection
3d1f02f9	294	TList collections_histoVtxBefSel;
	295	TList collections_histoVtxAftSel;
	296	TList collections_histoEtaBefSel;
	297	TList collections_histoEtaAftSel;
ae0fdd7d	298	TList collections_histoNChAftSel;
16ee3540	299	// TList collections_histoQVectorPos;
	300	// TList collections_histoQVectorNeg;
	301	TList collections_histoQVector;
	302	TList collections_histoEP;
c88234ad	303
	304	Int_t count = 0;
	305
	306	while ((obj = iter->Next())) {
	307	AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
	308	if (entry == 0)
	309	continue;
	310
	311	TH1I * histo = entry->GetHistoCuts();
	312	collections.Add(histo);
3d1f02f9	313	TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
	314	collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
	315	TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
	316	collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
	317	TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
	318	collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
	319	TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
	320	collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
ae0fdd7d	321	TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
ae0fdd7d	322	collections_histoNChAftSel.Add(histo_histoNChAftSel);
16ee3540	323	// TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
	324	// collections_histoQVectorPos.Add(histo_histoQVectorPos);
	325	// TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
	326	// collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
	327	TH1F * histo_histoQVector = entry->GetHistoQVector();
	328	collections_histoQVector.Add(histo_histoQVector);
	329	TH1F * histo_histoEP = entry->GetHistoEP();
	330	collections_histoEP.Add(histo_histoEP);
c88234ad	331	count++;
	332	}
	333
	334	fHistoCuts->Merge(&collections);
3d1f02f9	335	fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
	336	fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
	337	fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
	338	fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
ae0fdd7d	339	fHistoNChAftSel->Merge(&collections_histoNChAftSel);
16ee3540	340	// fHistoQVectorPos->Merge(&collections_histoQVectorPos);
	341	// fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
	342	fHistoQVector->Merge(&collections_histoQVector);
	343	fHistoEP->Merge(&collections_histoEP);
c88234ad	344
16ee3540	345
c88234ad	346	delete iter;
	347
	348	return count+1;
	349	}
	350