[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEextraCuts.cxx

/**************************************************************************
* Copyright(c) 1998-1999, 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.                  *
**************************************************************************/
/************************************************************************
 *                                                                      *
 * Extra cuts implemented by the ALICE Heavy Flavour Electron Group     *
 * Cuts stored here:                                                    *
 * - ITS pixels                                                         *
 * - TPC cluster ratio                                                  *
 * - TRD tracklets                                                      *
 *                                                                      *
 * Authors:                                                             *
 *   Markus Fasel <M.Fasel@gsi.de>                                      *
 *                                                                      *
 ************************************************************************/
#include <TAxis.h>
#include <TClass.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TString.h>

#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCParticle.h"

#include "AliHFEextraCuts.h"

ClassImp(AliHFEextraCuts)

//______________________________________________________
AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  AliCFCutBase(name, title),
  fCutCorrelation(0),
  fRequirements(0),
  fClusterRatioTPC(0.),
  fMinTrackletsTRD(0),
  fPixelITS(0),
  fQAlist(0x0),
  fDebugLevel(0)
{
  //
  // Default Constructor
  //
  memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
}

//______________________________________________________
AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  AliCFCutBase(c),
  fCutCorrelation(c.fCutCorrelation),
  fRequirements(c.fRequirements),
  fClusterRatioTPC(c.fClusterRatioTPC),
  fMinTrackletsTRD(c.fMinTrackletsTRD),
  fPixelITS(c.fPixelITS),
  fQAlist(0x0),
  fDebugLevel(c.fDebugLevel)
{
  //
  // Copy constructor
  // Performs a deep copy
  //
  memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
  if(IsQAOn()){
    fIsQAOn = kTRUE;
    fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
    fQAlist->SetOwner();
  }
}

//______________________________________________________
AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
  //
  // Assignment operator
  //
  if(this != &c){
    AliCFCutBase::operator=(c);
    fCutCorrelation = c.fCutCorrelation;
    fRequirements = c.fRequirements;
    fClusterRatioTPC = c.fClusterRatioTPC;
    fMinTrackletsTRD = c.fMinTrackletsTRD;
    fPixelITS = c.fPixelITS;
    fDebugLevel = c.fDebugLevel;

    memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
    if(IsQAOn()){
      fIsQAOn = kTRUE;
      fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
      fQAlist->SetOwner();
    }else fQAlist = 0x0;
  }
  return *this;
}

//______________________________________________________
AliHFEextraCuts::~AliHFEextraCuts(){
  //
  // Destructor
  //
  if(fQAlist){
    fQAlist->Delete();
    delete fQAlist;
  }
}

//______________________________________________________
Bool_t AliHFEextraCuts::IsSelected(TObject *o){
  //
  // Steering function for the track selection
  //
  if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
    return CheckESDCuts(dynamic_cast<AliESDtrack *>(o));
  }
  return CheckMCCuts(dynamic_cast<AliMCParticle *>(o));
}

//______________________________________________________
Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
  //
  // Checks cuts on reconstructed tracks
  // returns true if track is selected
  // QA histograms are filled before track selection and for
  // selected tracks after track selection
  //
  ULong64_t survivedCut = 0;	// Bitmap for cuts which are passed by the track, later to be compared with fRequirements
  if(IsQAOn()) FillQAhistosESD(track, kBeforeCuts);
  // Apply cuts
  Float_t impact_r, impact_z, ratioTPC;
  track->GetImpactParameters(impact_r, impact_z);
  // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
  ratioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
  UChar_t trdTracklets;
  #ifdef TRUNK
  trdTracklets = track->GetTRDntrackletsPID();
  #else
  trdTracklets = track->GetTRDpidQuality();
  #endif  
  UChar_t itsPixel = track->GetITSClusterMap();
  Int_t det, status1, status2;
  Float_t xloc, zloc;
  track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
  track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
  if(TESTBIT(fRequirements, kMinImpactParamR)){
    // cut on min. Impact Parameter in Radial direction
    if(impact_r >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
  }
  if(TESTBIT(fRequirements, kMinImpactParamZ)){
    // cut on min. Impact Parameter in Z direction
    if(impact_z >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
  }
  if(TESTBIT(fRequirements, kMaxImpactParamR)){
    // cut on max. Impact Parameter in Radial direction
    if(impact_r <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
  }
  if(TESTBIT(fRequirements, kMaxImpactParamZ)){
    // cut on max. Impact Parameter in Z direction
    if(impact_z <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
  }
  if(TESTBIT(fRequirements, kClusterRatioTPC)){
    // cut on min ratio of found TPC clusters vs findable TPC clusters
    if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
  }
  if(TESTBIT(fRequirements, kMinTrackletsTRD)){
    // cut on minimum number of TRD tracklets
    if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
  }
  if(TESTBIT(fRequirements, kPixelITS)){
    // cut on ITS pixel layers
    if(fDebugLevel > 0){
      printf("ITS cluster Map: ");
      PrintBitMap(itsPixel);
    }
    switch(fPixelITS){
      case kFirst: 
          if(itsPixel & BIT(0)) 
            SETBIT(survivedCut, kPixelITS);
          else if(!CheckITSstatus(status1))
            SETBIT(survivedCut, kPixelITS);
		    break;
      case kSecond: 
          if(itsPixel & BIT(1)) 
            SETBIT(survivedCut, kPixelITS);
          else if(!CheckITSstatus(status2))
            SETBIT(survivedCut, kPixelITS);
		    break;
      case kBoth: 
          if((itsPixel & BIT(0)) && (itsPixel & BIT(1))) 
            SETBIT(survivedCut, kPixelITS);
          else if(!CheckITSstatus(status1) && !CheckITSstatus(status2))
            SETBIT(survivedCut, kPixelITS);
		    break;
      case kAny: 
          if((itsPixel & BIT(0)) || (itsPixel & BIT(1))) 
            SETBIT(survivedCut, kPixelITS);
          else if(!CheckITSstatus(status1) || !CheckITSstatus(status2))
            SETBIT(survivedCut, kPixelITS);
		    break;
      default: break;
    }
    if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
  }
  if(fRequirements == survivedCut){
    //
    // Track selected
    //
    if(IsQAOn()) FillQAhistosESD(track, kAfterCuts);
    return kTRUE;
  }
  if(IsQAOn()) FillCutCorrelation(survivedCut);
  return kFALSE;
}

//______________________________________________________
Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle */*track*/){
  //
  // Checks cuts on Monte Carlo tracks
  // returns true if track is selected
  // QA histograms are filled before track selection and for
  // selected tracks after track selection
  //
  return kTRUE;	// not yet implemented
}

//______________________________________________________
void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
  //
  // Fill the QA histograms for ESD tracks
  // Function can be called before cuts or after cut application (second argument)
  //
  TList *container = dynamic_cast<TList *>(fQAlist->At(when));
  Float_t impact_r, impact_z;
  track->GetImpactParameters(impact_r, impact_z);
  (dynamic_cast<TH1F *>(container->At(0)))->Fill(impact_r);
  (dynamic_cast<TH1F *>(container->At(1)))->Fill(impact_z);
  // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
  (dynamic_cast<TH1F *>(container->At(2)))->Fill(track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.);
  #ifdef TRUNK
  (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
  #else
  (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDpidQuality());
  #endif
  UChar_t itsPixel = track->GetITSClusterMap();
  TH1 *pixelHist = dynamic_cast<TH1F *>(container->At(4));
  Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
  if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
    pixelHist->Fill(firstEntry + 3);
  else{
    if(itsPixel & BIT(0)){
      pixelHist->Fill(firstEntry);
      if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
      else pixelHist->Fill(firstEntry + 4);
    }
    if(itsPixel & BIT(1)){
      pixelHist->Fill(firstEntry + 1);
      if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
    }
  }
}

// //______________________________________________________
// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
//   //
//   // Fill the QA histograms for MC tracks
//   // Function can be called before cuts or after cut application (second argument)
//   // Not yet implemented
//   //
// }

//______________________________________________________
void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
  //
  // Fill cut correlation histograms for tracks that didn't pass cuts
  //
  TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2));
  for(Int_t icut = 0; icut < kNcuts; icut++){
    if(!TESTBIT(fRequirements, icut)) continue;
    for(Int_t jcut = icut; jcut < kNcuts; jcut++){
      if(!TESTBIT(fRequirements, jcut)) continue;
      if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
	correlation->Fill(icut, jcut);
    }
  }
}

//______________________________________________________
void AliHFEextraCuts::AddQAHistograms(TList *qaList){
  //
  // Add QA histograms
  // For each cut a histogram before and after track cut is created
  // Histos before respectively after cut are stored in different lists
  // Additionally a histogram with the cut correlation is created and stored
  // in the top directory
  //
  TList *histos[2];
  TH1 *histo1D = 0x0;
  TH2 *histo2D = 0x0;
  histos[0] = new TList();
  histos[0]->SetName("BeforeCut");
  histos[0]->SetOwner();
  histos[1] = new TList();
  histos[1]->SetName("AfterCut");
  histos[1]->SetOwner();
  TString cutstr[2] = {"before", "after"};
  for(Int_t icond = 0; icond < 2; icond++){
    histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamR%s", cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
    histo1D->GetXaxis()->SetTitle("Impact Parameter");
    histo1D->GetYaxis()->SetTitle("Number of Tracks");
    histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamZ%s", cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
    histo1D->GetXaxis()->SetTitle("Impact Parameter");
    histo1D->GetYaxis()->SetTitle("Number of Tracks");
    histos[icond]->AddAt((histo1D = new TH1F(Form("tpcClr%s", cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
    histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
    histo1D->GetYaxis()->SetTitle("Number of Tracks");
    histos[icond]->AddAt((histo1D = new TH1F(Form("trdTracklets%s", cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
    histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
    histo1D->GetYaxis()->SetTitle("Number of Tracks");
    histos[icond]->AddAt((histo1D = new TH1F(Form("itsPixel%s", cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
    histo1D->GetXaxis()->SetTitle("ITS Pixel");
    histo1D->GetYaxis()->SetTitle("Number of Tracks");
    Int_t first = histo1D->GetXaxis()->GetFirst();
    TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
    for(Int_t ilabel = 0; ilabel < 6; ilabel++)
      histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
  }
  fQAlist = new TList();
  fQAlist->SetOwner();
  fQAlist->SetName("HFelectronExtraCuts");
  fQAlist->AddAt(histos[0], 0);
  fQAlist->AddAt(histos[1], 1);
  // Add cut correlation
  fQAlist->AddAt((histo2D = new TH2F("cutcorrellation", "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
  TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
  Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
  for(Int_t icut = 0; icut < kNcuts; icut++){
    histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
    histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
  }
  qaList->AddLast(fQAlist);
}

//______________________________________________________
void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
  for(Int_t ibit = 32; ibit--; )
    printf("%d", bitmap & BIT(ibit) ? 1 : 0);
  printf("\n");
}

//______________________________________________________
Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus){
  //
  // Check whether ITS area is dead
  //
  Bool_t status;
  switch(itsStatus){
    case 2: status = kFALSE; break;
    case 3: status = kFALSE; break;
    case 7: status = kFALSE; break;
    default: status = kTRUE;
  }
  return status;
}
Commit	Line	Data
809a4336	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/************************************************************************
	16	* *
	17	* Extra cuts implemented by the ALICE Heavy Flavour Electron Group *
	18	* Cuts stored here: *
	19	* - ITS pixels *
	20	* - TPC cluster ratio *
	21	* - TRD tracklets *
	22	* *
	23	* Authors: *
	24	* Markus Fasel <M.Fasel@gsi.de> *
	25	* *
	26	************************************************************************/
	27	#include <TAxis.h>
	28	#include <TClass.h>
	29	#include <TH1F.h>
	30	#include <TH2F.h>
	31	#include <TList.h>
	32	#include <TString.h>
	33
	34	#include "AliESDtrack.h"
	35	#include "AliLog.h"
	36	#include "AliMCParticle.h"
	37
	38	#include "AliHFEextraCuts.h"
	39
	40	ClassImp(AliHFEextraCuts)
	41
	42	//______________________________________________________
	43	AliHFEextraCuts::AliHFEextraCuts(const Char_t name, const Char_t title):
	44	AliCFCutBase(name, title),
	45	fCutCorrelation(0),
	46	fRequirements(0),
	47	fClusterRatioTPC(0.),
	48	fMinTrackletsTRD(0),
	49	fPixelITS(0),
dbe3abbe	50	fQAlist(0x0),
dbe3abbe	51	fDebugLevel(0)
809a4336	52	{
	53	//
	54	// Default Constructor
	55	//
	56	memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
	57	}
	58
	59	//______________________________________________________
	60	AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
	61	AliCFCutBase(c),
	62	fCutCorrelation(c.fCutCorrelation),
	63	fRequirements(c.fRequirements),
	64	fClusterRatioTPC(c.fClusterRatioTPC),
	65	fMinTrackletsTRD(c.fMinTrackletsTRD),
	66	fPixelITS(c.fPixelITS),
dbe3abbe	67	fQAlist(0x0),
dbe3abbe	68	fDebugLevel(c.fDebugLevel)
809a4336	69	{
	70	//
	71	// Copy constructor
	72	// Performs a deep copy
	73	//
	74	memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
	75	if(IsQAOn()){
	76	fIsQAOn = kTRUE;
	77	fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
	78	fQAlist->SetOwner();
	79	}
	80	}
	81
	82	//______________________________________________________
	83	AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
	84	//
	85	// Assignment operator
	86	//
	87	if(this != &c){
	88	AliCFCutBase::operator=(c);
	89	fCutCorrelation = c.fCutCorrelation;
	90	fRequirements = c.fRequirements;
	91	fClusterRatioTPC = c.fClusterRatioTPC;
	92	fMinTrackletsTRD = c.fMinTrackletsTRD;
	93	fPixelITS = c.fPixelITS;
dbe3abbe	94	fDebugLevel = c.fDebugLevel;
809a4336	95
	96	memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
	97	if(IsQAOn()){
	98	fIsQAOn = kTRUE;
	99	fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
	100	fQAlist->SetOwner();
	101	}else fQAlist = 0x0;
	102	}
	103	return *this;
	104	}
	105
	106	//______________________________________________________
	107	AliHFEextraCuts::~AliHFEextraCuts(){
	108	//
	109	// Destructor
	110	//
	111	if(fQAlist){
	112	fQAlist->Delete();
	113	delete fQAlist;
	114	}
	115	}
	116
	117	//______________________________________________________
	118	Bool_t AliHFEextraCuts::IsSelected(TObject *o){
	119	//
	120	// Steering function for the track selection
	121	//
	122	if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
	123	return CheckESDCuts(dynamic_cast<AliESDtrack *>(o));
	124	}
	125	return CheckMCCuts(dynamic_cast<AliMCParticle *>(o));
	126	}
	127
	128	//______________________________________________________
	129	Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
	130	//
	131	// Checks cuts on reconstructed tracks
	132	// returns true if track is selected
	133	// QA histograms are filled before track selection and for
	134	// selected tracks after track selection
	135	//
	136	ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
	137	if(IsQAOn()) FillQAhistosESD(track, kBeforeCuts);
	138	// Apply cuts
	139	Float_t impact_r, impact_z, ratioTPC;
	140	track->GetImpactParameters(impact_r, impact_z);
	141	// printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
	142	ratioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
	143	UChar_t trdTracklets;
	144	#ifdef TRUNK
	145	trdTracklets = track->GetTRDntrackletsPID();
	146	#else
	147	trdTracklets = track->GetTRDpidQuality();
	148	#endif
	149	UChar_t itsPixel = track->GetITSClusterMap();
dbe3abbe	150	Int_t det, status1, status2;
	151	Float_t xloc, zloc;
	152	track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
	153	track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
809a4336	154	if(TESTBIT(fRequirements, kMinImpactParamR)){
	155	// cut on min. Impact Parameter in Radial direction
	156	if(impact_r >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
	157	}
	158	if(TESTBIT(fRequirements, kMinImpactParamZ)){
	159	// cut on min. Impact Parameter in Z direction
	160	if(impact_z >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
	161	}
	162	if(TESTBIT(fRequirements, kMaxImpactParamR)){
	163	// cut on max. Impact Parameter in Radial direction
	164	if(impact_r <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
	165	}
	166	if(TESTBIT(fRequirements, kMaxImpactParamZ)){
	167	// cut on max. Impact Parameter in Z direction
	168	if(impact_z <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
	169	}
	170	if(TESTBIT(fRequirements, kClusterRatioTPC)){
	171	// cut on min ratio of found TPC clusters vs findable TPC clusters
	172	if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
	173	}
	174	if(TESTBIT(fRequirements, kMinTrackletsTRD)){
	175	// cut on minimum number of TRD tracklets
	176	if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
	177	}
	178	if(TESTBIT(fRequirements, kPixelITS)){
	179	// cut on ITS pixel layers
dbe3abbe	180	if(fDebugLevel > 0){
	181	printf("ITS cluster Map: ");
	182	PrintBitMap(itsPixel);
	183	}
809a4336	184	switch(fPixelITS){
dbe3abbe	185	case kFirst:
	186	if(itsPixel & BIT(0))
	187	SETBIT(survivedCut, kPixelITS);
	188	else if(!CheckITSstatus(status1))
	189	SETBIT(survivedCut, kPixelITS);
	190	break;
	191	case kSecond:
	192	if(itsPixel & BIT(1))
	193	SETBIT(survivedCut, kPixelITS);
	194	else if(!CheckITSstatus(status2))
	195	SETBIT(survivedCut, kPixelITS);
	196	break;
	197	case kBoth:
	198	if((itsPixel & BIT(0)) && (itsPixel & BIT(1)))
	199	SETBIT(survivedCut, kPixelITS);
	200	else if(!CheckITSstatus(status1) && !CheckITSstatus(status2))
	201	SETBIT(survivedCut, kPixelITS);
	202	break;
	203	case kAny:
	204	if((itsPixel & BIT(0)) \|\| (itsPixel & BIT(1)))
	205	SETBIT(survivedCut, kPixelITS);
	206	else if(!CheckITSstatus(status1) \|\| !CheckITSstatus(status2))
	207	SETBIT(survivedCut, kPixelITS);
	208	break;
809a4336	209	default: break;
809a4336	210	}
dbe3abbe	211	if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
809a4336	212	}
	213	if(fRequirements == survivedCut){
	214	//
	215	// Track selected
	216	//
	217	if(IsQAOn()) FillQAhistosESD(track, kAfterCuts);
	218	return kTRUE;
	219	}
	220	if(IsQAOn()) FillCutCorrelation(survivedCut);
	221	return kFALSE;
	222	}
	223
	224	//______________________________________________________
	225	Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle /track*/){
	226	//
	227	// Checks cuts on Monte Carlo tracks
	228	// returns true if track is selected
	229	// QA histograms are filled before track selection and for
	230	// selected tracks after track selection
	231	//
	232	return kTRUE; // not yet implemented
	233	}
	234
	235	//______________________________________________________
	236	void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
	237	//
	238	// Fill the QA histograms for ESD tracks
	239	// Function can be called before cuts or after cut application (second argument)
	240	//
	241	TList container = dynamic_cast<TList >(fQAlist->At(when));
	242	Float_t impact_r, impact_z;
	243	track->GetImpactParameters(impact_r, impact_z);
	244	(dynamic_cast<TH1F *>(container->At(0)))->Fill(impact_r);
	245	(dynamic_cast<TH1F *>(container->At(1)))->Fill(impact_z);
	246	// printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
	247	(dynamic_cast<TH1F *>(container->At(2)))->Fill(track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.);
	248	#ifdef TRUNK
	249	(dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
	250	#else
	251	(dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDpidQuality());
	252	#endif
	253	UChar_t itsPixel = track->GetITSClusterMap();
	254	TH1 pixelHist = dynamic_cast<TH1F >(container->At(4));
	255	Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
	256	if(!((itsPixel & BIT(0)) \|\| (itsPixel & BIT(1))))
	257	pixelHist->Fill(firstEntry + 3);
	258	else{
	259	if(itsPixel & BIT(0)){
	260	pixelHist->Fill(firstEntry);
	261	if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
	262	else pixelHist->Fill(firstEntry + 4);
	263	}
	264	if(itsPixel & BIT(1)){
	265	pixelHist->Fill(firstEntry + 1);
	266	if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
	267	}
	268	}
	269	}
	270
	271	// //______________________________________________________
	272	// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
	273	// //
	274	// // Fill the QA histograms for MC tracks
	275	// // Function can be called before cuts or after cut application (second argument)
276	// // Not yet implemented
277	// //
278	// }
279
280	//______________________________________________________
281	void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
282	//
283	// Fill cut correlation histograms for tracks that didn't pass cuts
284	//
285	TH2 correlation = dynamic_cast<TH2F >(fQAlist->At(2));
286	for(Int_t icut = 0; icut < kNcuts; icut++){
287	if(!TESTBIT(fRequirements, icut)) continue;
288	for(Int_t jcut = icut; jcut < kNcuts; jcut++){
289	if(!TESTBIT(fRequirements, jcut)) continue;
290	if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
291	correlation->Fill(icut, jcut);
292	}
293	}
294	}
295
296	//______________________________________________________
297	void AliHFEextraCuts::AddQAHistograms(TList *qaList){
298	//
299	// Add QA histograms
300	// For each cut a histogram before and after track cut is created
301	// Histos before respectively after cut are stored in different lists
302	// Additionally a histogram with the cut correlation is created and stored
303	// in the top directory
304	//
305	TList *histos[2];
306	TH1 *histo1D = 0x0;
307	TH2 *histo2D = 0x0;
308	histos[0] = new TList();
309	histos[0]->SetName("BeforeCut");
310	histos[0]->SetOwner();
311	histos[1] = new TList();
312	histos[1]->SetName("AfterCut");
313	histos[1]->SetOwner();
314	TString cutstr[2] = {"before", "after"};
315	for(Int_t icond = 0; icond < 2; icond++){
316	histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamR%s", cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
317	histo1D->GetXaxis()->SetTitle("Impact Parameter");
318	histo1D->GetYaxis()->SetTitle("Number of Tracks");
319	histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamZ%s", cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
320	histo1D->GetXaxis()->SetTitle("Impact Parameter");
321	histo1D->GetYaxis()->SetTitle("Number of Tracks");
322	histos[icond]->AddAt((histo1D = new TH1F(Form("tpcClr%s", cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
323	histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
324	histo1D->GetYaxis()->SetTitle("Number of Tracks");
325	histos[icond]->AddAt((histo1D = new TH1F(Form("trdTracklets%s", cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
326	histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
327	histo1D->GetYaxis()->SetTitle("Number of Tracks");
328	histos[icond]->AddAt((histo1D = new TH1F(Form("itsPixel%s", cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
329	histo1D->GetXaxis()->SetTitle("ITS Pixel");
330	histo1D->GetYaxis()->SetTitle("Number of Tracks");
331	Int_t first = histo1D->GetXaxis()->GetFirst();
332	TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
333	for(Int_t ilabel = 0; ilabel < 6; ilabel++)
334	histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
335	}
336	fQAlist = new TList();
337	fQAlist->SetOwner();
338	fQAlist->SetName("HFelectronExtraCuts");
339	fQAlist->AddAt(histos[0], 0);
340	fQAlist->AddAt(histos[1], 1);
341	// Add cut correlation
342	fQAlist->AddAt((histo2D = new TH2F("cutcorrellation", "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
343	TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
344	Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
345	for(Int_t icut = 0; icut < kNcuts; icut++){
346	histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
347	histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
348	}
349	qaList->AddLast(fQAlist);
350	}
dbe3abbe	351
	352	//______________________________________________________
	353	void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
	354	for(Int_t ibit = 32; ibit--; )
	355	printf("%d", bitmap & BIT(ibit) ? 1 : 0);
	356	printf("\n");
	357	}
	358
	359	//______________________________________________________
	360	Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus){
	361	//
	362	// Check whether ITS area is dead
	363	//
	364	Bool_t status;
	365	switch(itsStatus){
	366	case 2: status = kFALSE; break;
	367	case 3: status = kFALSE; break;
	368	case 7: status = kFALSE; break;
	369	default: status = kTRUE;
	370	}
	371	return status;
	372	}