[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 "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),fIsMC(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorPosCutMin(0), fQVectorPosCutMax(0), fQVectorNegCutMin(0), fQVectorNegCutMax(0), fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVectorPos(0),fHistoQVectorNeg(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",3000,-0.5,2999.5);
  fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
  fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,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;
}

//______________________________________________________
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
  }
  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 (!fTrackCuts->IsSelected(aodTrack,kFALSE)) 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;
  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;
}
//______________________________________________________

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;

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