[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), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(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;
  fVertexCutMin=-10.0;
  fVertexCutMax=10.0;
  fMultiplicityCutMin=0.0;
  fMultiplicityCutMax=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() && CheckMultiplicityCut()){ //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++;
    }
  }
  //Printf("NCHARGED_EvSel : %d",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 (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
    {
      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::CheckMultiplicityCut()
{
  // Check multiplicity cut
  Int_t Ncharged=0;
  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
    AliAODTrack* track = fAOD->GetTrack(iTracks);
    if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
    Ncharged++;
  }
  //Printf("NCHARGED_cut : %d",Ncharged);
  if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
  
  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 << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
  cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< 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
10a8ccbe	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), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(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;
10a8ccbe	67	fVertexCutMin=-10.0;
10a8ccbe	68	fVertexCutMax=10.0;
351d117d	69	fMultiplicityCutMin=0.0;
10a8ccbe	70	fMultiplicityCutMax=10000.0;
829b5a81	71	fTrackBits=1;
c88234ad	72	}
	73
	74	//______________________________________________________
f6a38178	75	Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts *trackcuts)
c88234ad	76	{
f6a38178	77	// Returns true if Event Cuts are selected and applied
	78	fAOD = aod;
	79	fTrackCuts = trackcuts;
	80	fHistoCuts->Fill(kProcessedEvents);
decf69d9	81	Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
decf69d9	82	if(!IsPhysSel)return IsPhysSel;
8961bb9a	83	fHistoCuts->Fill(kPhysSelEvents);
f6a38178	84	//loop on tracks, before event selection, filling QA histos
	85	AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
	86	if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
decf69d9	87	fIsSelected =kFALSE;
10a8ccbe	88	if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
00493191	89	fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
decf69d9	90	}
f6a38178	91	if(fIsSelected){
	92	fHistoCuts->Fill(kAcceptedEvents);
	93	if(vertex)fHistoVtxAftSel->Fill(vertex->GetZ());
	94	}
ae0fdd7d	95	Int_t Nch=0;
f6a38178	96	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
f6a38178	97	AliAODTrack* track = fAOD->GetTrack(iTracks);
decf69d9	98	if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
f6a38178	99	fHistoEtaBefSel->Fill(track->Eta());
ae0fdd7d	100	if(fIsSelected){
	101	fHistoEtaAftSel->Fill(track->Eta());
	102	Nch++;
	103	}
f6a38178	104	}
194d6dce	105	//Printf("NCHARGED_EvSel : %d",Nch);
ae0fdd7d	106	if(fIsSelected)fHistoNChAftSel->Fill(Nch);
f6a38178	107	return fIsSelected;
c88234ad	108	}
	109
	110	//______________________________________________________
	111	Bool_t AliSpectraAODEventCuts::CheckVtxRange()
	112	{
	113	// reject events outside of range
f6a38178	114	AliAODVertex * vertex = fAOD->GetPrimaryVertex();
	115	if (!vertex)
	116	{
c88234ad	117	fHistoCuts->Fill(kVtxNoEvent);
c88234ad	118	return kFALSE;
f6a38178	119	}
10a8ccbe	120	if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
f6a38178	121	{
c88234ad	122	return kTRUE;
f6a38178	123	}
	124	fHistoCuts->Fill(kVtxRange);
	125	return kFALSE;
c88234ad	126	}
	127
	128	//______________________________________________________
	129	Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
	130	{
f6a38178	131	// Check centrality cut
972a21ad	132	Double_t cent=0;
936e2842	133	if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile("V0M");
972a21ad	134	else cent=ApplyCentralityPatchAOD049();
	135
	136	if ( (cent <= fCentralityCutMax) && (cent >= fCentralityCutMin) ) return kTRUE;
f6a38178	137	fHistoCuts->Fill(kVtxCentral);
f6a38178	138	return kFALSE;
c88234ad	139	}
c88234ad	140
10a8ccbe	141	//______________________________________________________
	142	Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
	143	{
	144	// Check multiplicity cut
	145	Int_t Ncharged=0;
	146	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
	147	AliAODTrack* track = fAOD->GetTrack(iTracks);
	148	if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
	149	Ncharged++;
	150	}
194d6dce	151	//Printf("NCHARGED_cut : %d",Ncharged);
10a8ccbe	152	if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
	153
	154	return kFALSE;
	155	}
	156
decf69d9	157	//______________________________________________________
	158	Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
	159	{
	160	// Check qvector cut
	161	/// FIXME: Q vector
	162	// //Selection on QVector, before ANY other selection on the event
	163	// //Spectra MUST be normalized wrt events AFTER the selection on Qvector
16ee3540	164	// Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
	165	// Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
	166
	167	// Int_t multPos = 0;
	168	// Int_t multNeg = 0;
	169	// for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
	170	// AliAODTrack* aodTrack = fAOD->GetTrack(iT);
	171	// if (!aodTrack->TestFilterBit(fTrackBits)) continue;
	172	// if (aodTrack->Eta() >= 0){
	173	// multPos++;
	174	// Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
	175	// Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
	176	// } else {
	177	// multNeg++;
	178	// Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
	179	// Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
	180	// }
	181	// }
	182	// Double_t qPos=-999;
	183	// if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPosQx2EtaPos + Qy2EtaPosQy2EtaPos)/multPos);
	184	// Double_t qNeg=-999;
	185	// if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNegQx2EtaNeg + Qy2EtaNegQy2EtaNeg)/multNeg);
	186	//if(qPos<fQVectorPosCutMin \|\| qPos>fQVectorPosCutMax \|\| qNeg<fQVectorNegCutMin \|\| qNeg>fQVectorNegCutMax)return kFALSE;
decf69d9	187
16ee3540	188	Double_t qxEPVZERO = 0, qyEPVZERO = 0;
	189	Double_t qVZERO = -999;
	190	Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
00493191	191
16ee3540	192	qVZERO= TMath::Sqrt(qxEPVZEROqxEPVZERO + qyEPVZEROqyEPVZERO);
	193	if(qVZERO<fQVectorCutMin \|\| qVZERO>fQVectorCutMax)return kFALSE;
	194	fHistoQVector->Fill(qVZERO);
	195	fHistoEP->Fill(psi);
00493191	196
decf69d9	197	fHistoCuts->Fill(kQVector);
16ee3540	198	// fHistoQVectorPos->Fill(qPos);
16ee3540	199	// fHistoQVectorNeg->Fill(qNeg);
decf69d9	200	return kTRUE;
decf69d9	201	}
decf69d9	202
c88234ad	203	//______________________________________________________
	204	void AliSpectraAODEventCuts::PrintCuts()
	205	{
f6a38178	206	// print info about event cuts
	207	cout << "Event Stats" << endl;
	208	cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
	209	cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
936e2842	210	cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
f6a38178	211	cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
	212	cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
	213	cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
decf69d9	214	cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
829b5a81	215	cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
16ee3540	216	// cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
	217	// cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
	218	cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
10a8ccbe	219	cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
10a8ccbe	220	cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
829b5a81	221	cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
f6a38178	222	}
c88234ad	223	//______________________________________________________
c88234ad	224
972a21ad	225	Double_t AliSpectraAODEventCuts::ApplyCentralityPatchAOD049()
	226	{
	227	//
	228	//Apply centrality patch for AOD049 outliers
	229	//
	230	// if (fCentralityType!="V0M") {
	231	// AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
	232	// return -999.0;
	233	// }
	234	AliCentrality *centrality = fAOD->GetCentrality();
	235	if (!centrality) {
	236	AliWarning("Cannot get centrality from AOD event.");
	237	return -999.0;
	238	}
	239
	240	Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
	241	/*
	242	Bool_t isSelRun = kFALSE;
	243	Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
	244	if(cent<0){
	245	Int_t quality = centrality->GetQuality();
	246	if(quality<=1){
	247	cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
	248	} else {
	249	Int_t runnum=aodEvent->GetRunNumber();
	250	for(Int_t ir=0;ir<5;ir++){
	251	if(runnum==selRun[ir]){
	252	isSelRun=kTRUE;
	253	break;
	254	}
	255	}
	256	if((quality==8\|\|quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
	257	}
	258	}
	259	*/
	260	if(cent>=0.0) {
	261	Float_t v0 = 0.0;
	262	AliAODEvent aodEvent = (AliAODEvent )fAOD;
	263	AliAODVZERO aodV0 = (AliAODVZERO ) aodEvent->GetVZEROData();
	264	v0+=aodV0->GetMTotV0A();
	265	v0+=aodV0->GetMTotV0C();
	266	if ( (cent==0) && (v0<19500) ) {
	267	AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
	268	return -999.0;
	269	}
	270	Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
	271	Float_t val = 1.30552 + 0.147931 * v0;
	272
	273	Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
	274	120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
	275	92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
	276	68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
	277	51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
	278	37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
	279	26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
	280	19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
	281	12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
	282	13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
	283	};
	284
	285	if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
	286	AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
	287	return -999.0;
	288	}
289	} else {
290	//force it to be -999. whatever the negative value was
291	cent = -999.;
292	}
293	return cent;
294	}
295
296	//______________________________________________________
297
298
c88234ad	299	Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
	300	{
	301	// Merge a list of AliSpectraAODEventCuts objects with this.
	302	// Returns the number of merged objects (including this).
	303
	304	// AliInfo("Merging");
	305
	306	if (!list)
	307	return 0;
	308
	309	if (list->IsEmpty())
	310	return 1;
	311
	312	TIterator* iter = list->MakeIterator();
	313	TObject* obj;
	314
	315	// collections of all histograms
ae0fdd7d	316	TList collections;//FIXME we should use only 1 collection
3d1f02f9	317	TList collections_histoVtxBefSel;
	318	TList collections_histoVtxAftSel;
	319	TList collections_histoEtaBefSel;
	320	TList collections_histoEtaAftSel;
ae0fdd7d	321	TList collections_histoNChAftSel;
16ee3540	322	// TList collections_histoQVectorPos;
	323	// TList collections_histoQVectorNeg;
	324	TList collections_histoQVector;
	325	TList collections_histoEP;
c88234ad	326
	327	Int_t count = 0;
	328
	329	while ((obj = iter->Next())) {
	330	AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
	331	if (entry == 0)
	332	continue;
	333
	334	TH1I * histo = entry->GetHistoCuts();
	335	collections.Add(histo);
3d1f02f9	336	TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
	337	collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
	338	TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
	339	collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
	340	TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
	341	collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
	342	TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
	343	collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
ae0fdd7d	344	TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
ae0fdd7d	345	collections_histoNChAftSel.Add(histo_histoNChAftSel);
16ee3540	346	// TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
	347	// collections_histoQVectorPos.Add(histo_histoQVectorPos);
	348	// TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
	349	// collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
	350	TH1F * histo_histoQVector = entry->GetHistoQVector();
	351	collections_histoQVector.Add(histo_histoQVector);
	352	TH1F * histo_histoEP = entry->GetHistoEP();
	353	collections_histoEP.Add(histo_histoEP);
c88234ad	354	count++;
	355	}
	356
	357	fHistoCuts->Merge(&collections);
3d1f02f9	358	fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
	359	fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
	360	fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
	361	fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
ae0fdd7d	362	fHistoNChAftSel->Merge(&collections_histoNChAftSel);
16ee3540	363	// fHistoQVectorPos->Merge(&collections_histoQVectorPos);
	364	// fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
	365	fHistoQVector->Merge(&collections_histoQVector);
	366	fHistoEP->Merge(&collections_histoEP);
c88234ad	367
16ee3540	368
c88234ad	369	delete iter;
	370
	371	return count+1;
	372	}
	373