[u/mrichter/AliRoot.git] / PWGUD / DIFFRACTIVE / xsAndTwoProng / AliAnalysisTaskCDskimESD.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.                  *
**************************************************************************/
//
// $Id$
//
// Task to skim ESD files, basic idea found at PWGGA/EMCALTasks as
// AliEsdSkimTask. This version is modified to produce a exact copy of the ESD
// stream containing only events which are of interest for the central
// diffractive analysis (PWGUD) and a control sample. In order to be able to
// process MC data as well, the usual concept using an output slot is for
// output data is violated and an output file is directly created without
// a corresponding output slot.
//
// Author:
//  Felix Reidt <Felix.Reidt@cern.ch>

// ROOT headers
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TAxis.h>
#include <TRandom3.h>
//#include <TKey.h>
//#include <TROOT.h>

// AliRoot headers
#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliESDFMD.h"
#include "AliEventplane.h"
#include "AliInputEventHandler.h"
#include "AliMultiplicity.h"
#include "AliVEventHandler.h"
#include "AliMCEventHandler.h"
//#include "AliMCEvent.h"
//#include "AliStack.h"
//#include "AliGenEventHeader.h"
//#include "AliRunLoader.h"

// own classes
#include "AliCDMesonBase.h"
#include "AliCDMesonUtils.h"

// own header
#include "AliAnalysisTaskCDskimESD.h"

ClassImp(AliAnalysisTaskCDskimESD)

//______________________________________________________________________________
AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD(const char *opt,
                                                   Bool_t reduceGapEvents)
	: AliAnalysisTaskSE(opt)
	, fInputTree(0x0)
	, fInputTreeCopy(0x0)
	, fOutputTree(0x0)
	, fDoMC(kFALSE)
	, fWrkDir(0x0)
	, fKinematicsFile(0x0)
	, fTrackRefsFile(0x0)
	, fTEinput(0x0)
	, fTE(0x0)
	, fgaliceFile(0x0)
	, fDoGapCond(kTRUE)
	, fCurrentGapCondition(0)
	, fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
	                        // only events containing a V0 double gap or
	, fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
	, fReduceGapEvents(reduceGapEvents)
	, fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
	                       AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
	, fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
	, fStatsFlow(0x0)
	, fSkimmingList(0x0)
	, fFileName()
	, fEventNumberInFile(-1)
	, fRunNumber(-1)
	, fEventTime(0)
{
	//
	// Constructor.
	//

	if (!opt)
		return;

	//fBranchNames = "ESD:AliESDHeader.,AliESDRun."; // TODO don't we need all?

	// check whether we are running on MC data
	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	AliMCEventHandler* mcH =
		dynamic_cast<AliMCEventHandler*>(am->GetMCtruthEventHandler());
	if(mcH) {
		fDoMC = kTRUE;

		fWrkDir = gDirectory; // save the old directory
		fKinematicsFile = new TFile("Kinematics.root", "recreate");
		fKinematicsFile->Close();
		delete fKinematicsFile;
		fKinematicsFile = 0x0;
		gDirectory = fWrkDir; // restore the old directory

		fTrackRefsFile = new TFile("TrackRefs.root", "recreate");
		fTrackRefsFile->Close();
		delete fTrackRefsFile;
		fTrackRefsFile = 0x0;
		gDirectory = fWrkDir; // restore the old directory
	}

	DefineOutput(1, TTree::Class());
	DefineOutput(2, TH1::Class());
	DefineOutput(3, TTree::Class());
	if (fDoMC) DefineOutput(4, TTree::Class());
}


//______________________________________________________________________________
AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD()
	: AliAnalysisTaskSE()
	, fInputTree(0x0)
	, fInputTreeCopy(0x0)
	, fOutputTree(0x0)
	, fDoMC(kFALSE)
	, fWrkDir(0x0)
	, fKinematicsFile(0x0)
	, fTrackRefsFile(0x0)
	, fTEinput(0x0)
	, fTE(0x0)
	, fgaliceFile(0x0)
	, fDoGapCond(kTRUE)
	, fCurrentGapCondition(0)
	, fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
	                        // only events containing a V0 double gap or
	, fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
	, fReduceGapEvents(kFALSE)
	, fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
	                       AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
	, fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
	, fStatsFlow(0x0)
	, fSkimmingList(0x0)
	, fFileName()
	, fEventNumberInFile(-1)
	, fRunNumber(-1)
	, fEventTime(0)
{
	//
	// Default Constructor
	//
}


//______________________________________________________________________________
AliAnalysisTaskCDskimESD::~AliAnalysisTaskCDskimESD()
{
	//
	// deconstructor
	//

	/* // this delete lead to chrases caused by the garbage collection ...
	if (fOutputTree) {
		delete fOutputTree;
		fOutputTree = 0x0;
		}
	if (fInputTreeCopy) {
		delete fInputTreeCopy;
		fInputTreeCopy = 0x0;
	}
	if (fStatsFlow) {
		delete fStatsFlow;
		fStatsFlow = 0x0;
	} */
	if (fDoMC && fKinematicsFile) {
		fKinematicsFile->Close();
		delete fKinematicsFile;
		fKinematicsFile = 0x0;
	}
	if (fDoMC && fTrackRefsFile) {
		fTrackRefsFile->Close();
		delete fTrackRefsFile;
		fTrackRefsFile = 0x0;
	}
}


//______________________________________________________________________________
void AliAnalysisTaskCDskimESD::UserExec(Option_t */*opt*/)
{
	//
  // Process event.
	//

	fStatsFlow->Fill(0.); // UserExec(...) <= bin name

	if (fDoMC) {
		// replace dummy tree by the correct tree
		AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
		if(!fTEinput) {
			fTEinput =
				((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
			if (fTEinput) {
				delete fTE;
				fTE = fTEinput->CloneTree(0);
				fTE->SetDirectory(fgaliceFile);
				fTE->SetAutoFlush(-10*1024*1024);
				//CopygALICE(am);
			}
		}
		else if (fTEinput != am->GetMCtruthEventHandler()->GetTree()) {
			fTEinput =
				((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
			fTEinput->CopyAddresses(fTE);
		}
	}

	// check whether the ESD input tree is still the same - NOT NECESSARY
	//TTree* currentInput =
	//	AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree();
	//if (fInputTree != currentInput) { // ESD input tree has changed
	//	puts("CHANGED INPUT TREE");
	//	fInputTree = currentInput;
	//	fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");
	//	fInputTreeCopy->CopyAddresses(fOutputTree);
	//}

	AliESDEvent *esdin = dynamic_cast<AliESDEvent*>(InputEvent());
	if (!esdin || !fInputTree || !fInputTreeCopy || !fOutputTree ||
	    (fDoMC && (!fTEinput || !fTE))) {
		// check whether everything is ready
		PostData(1, fOutputTree);
		PostData(2, fStatsFlow);
		PostData(3, fSkimmingList);
		if (fDoMC) PostData(4, fTE);
		return;
	}

	fStatsFlow->Fill(1.); // Ready <= bin name

	//============================================================================
	// event selection, every event passing the next lines is written to the new
	// ESD stream
	//
	// modify this lines in order to adjust the ESD selection to your needs

	// determine current gap condition
	fCurrentGapCondition = AliCDMesonUtils::GetGapConfig(esdin, 0x0, 0x0, 0x0,
	                                                     0x0, 0x0, 0x0, 0x0, 0x0);

	TRandom3 rnd(0);
	if (!(AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) ||
	     (fRequestedGapCondition & fCurrentGapCondition) || fReduceGapEvents) {
		// check whether the event satifies the required gap condition and whether
		// all double-gap events will be stored

		if (fReduceGapEvents &&
		    (AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) &&
		    !(fRequestedGapCondition & fCurrentGapCondition)) {
			// not all double gap events can be stored, but the gap determination
			// was successful and at least a according to the less stringent condition
			// a double gap was found
			if (!(fRefinedGapCondition & fCurrentGapCondition)) {
				// event fulfilled the more stringent gap and hence is will be stored
				fStatsFlow->Fill(5.);
			}
			else if (rnd.Rndm() > (1.-fReducedGapFraction)) {
				// randomly select a fraction of
				fStatsFlow->Fill(4.); // Double Gap Event <= bin name
			}
			else {
				// reject event (no refined gap and not randomly selected)
				fStatsFlow->Fill(3.); // Event Rejected <= bin name
				PostData(1, fOutputTree);
				PostData(2, fStatsFlow);
				PostData(3, fSkimmingList);
				if (fDoMC) PostData(4, fTE);
				return;
			}
		}
		else if (rnd.Rndm() > (1.-fNoGapFraction)) { // randomly selected
			fStatsFlow->Fill(2.); // Control Sample Event <= bin name
		}
		else {
			fStatsFlow->Fill(3.); // Event Rejected <= bin name
			PostData(1, fOutputTree);
			PostData(2, fStatsFlow);
			PostData(3, fSkimmingList);
			if (fDoMC) PostData(4, fTE);
			return;
		}
	}
	else {
		fStatsFlow->Fill(4.); // Double Gap Event <= bin name
	}
	// end of event selection
	//============================================================================

	// load the current event in the cloned input tree
	//
	// unfortunately the Entry() gives the entry number with in the current input
	// tree, not within the chain for MC data
	// => use GetReadEntry() of the ESD tree (value return by Entry() for real
	//data)
	fInputTreeCopy->GetEvent(fInputTree->GetReadEntry());
	fOutputTree->Fill(); // fill the current event into the new ESD stream
	//printf("Entry()=%lld\nEntries=%lld\nChainOffset=%lld\nReadEntry=%lld\n",
	//       Entry(), fInputTreeCopy->GetEntries(),
	//       fInputTreeCopy->GetChainOffset(), fInputTree->GetReadEntry());

	// MC specific stuff ---------------------------------------------------------
	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

	if (fDoMC) { // do the MC related part of the skimming
		Long64_t iEntry = fOutputTree->GetEntries() - 1;
		AliMCEventHandler* mcHandler =
			(AliMCEventHandler*)am->GetMCtruthEventHandler();
		TTree* treeK = mcHandler->TreeK();
		TTree* treeTR = mcHandler->TreeTR();
		if (!treeK) {
			AliFatal("TreeK not found!");
			return;
		}
		if (!treeTR) {
			AliFatal("TreeTR not found!");
			return;
		}

		fWrkDir = gDirectory;
		gDirectory = fKinematicsFile;
		TString dirName = TString(Form("Event%ld", (Long_t)iEntry));
		gDirectory->mkdir(dirName);
		gDirectory->cd(dirName);
		TTree* outputTreeK = treeK->CloneTree(); // copy the whole tree
		if (!outputTreeK) {
			AliFatal("Could not clone TreeK!");
			return;
		}
		outputTreeK->Write();
		treeK->CopyAddresses(outputTreeK, kTRUE); // separate clone again
		treeK->GetListOfClones()->Remove((TObject*)outputTreeK);
		outputTreeK->Delete();
		outputTreeK = 0x0;

		gDirectory = fTrackRefsFile;
		gDirectory->mkdir(dirName);
		gDirectory->cd(dirName);
		TTree* outputTreeTR = treeTR->CloneTree(); // copy the whole tree
		if (!outputTreeTR) {
			AliFatal("Could not clone TreeTR");
			return;
		}
		outputTreeTR->Write();
		treeTR->CopyAddresses(outputTreeTR, kTRUE); // separate clone again
		treeTR->GetListOfClones()->Remove((TObject*)outputTreeTR);
		outputTreeTR->Delete();
		outputTreeTR = 0x0;
		gDirectory = fWrkDir;

		TTree* currentTEinput =
			((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
		if (fTEinput != currentTEinput) {
			// when the input file is changed, the kinematics tree address changes
			// hence our output tree fTE has to be linked to the new one!
			fTEinput = currentTEinput;
			fTEinput->CopyAddresses(fTE);
		}
		fTE->Fill(); // fill MC event headers
	}

	fStatsFlow->Fill(6.); // Tree Filled <= bin name

	// remember which event we stored
	fEventNumberInFile = esdin->GetEventNumberInFile();
	fRunNumber = esdin->GetRunNumber();
	fFileName = fInputTree->GetCurrentFile()->GetName();
	fEventTime = esdin->GetTimeStamp();
	fSkimmingList->Fill();

	PostData(1, fOutputTree);
	PostData(2, fStatsFlow);
	PostData(3, fSkimmingList);
	if (fDoMC) PostData(4, fTE);
}

//______________________________________________________________________________
void AliAnalysisTaskCDskimESD::UserCreateOutputObjects()
{
	// Create output objects.
	fStatsFlow = (TH1*)(new TH1F("skimmingStatsFlow", "", 7, 0., 7.));
	TAxis* x = fStatsFlow->GetXaxis();
	x->SetBinLabel(1, "UserExec(...)");
	x->SetBinLabel(2, "Ready");
	x->SetBinLabel(3, "Control Sample Event");
	x->SetBinLabel(4, "Rejected Event");
	x->SetBinLabel(5, "Double Gap Event");
	x->SetBinLabel(6, "Refined Double Gap Event");
	x->SetBinLabel(7, "Tree Filled");

	PostData(2, fStatsFlow);

	// TODO implement fSkimmingList
	fSkimmingList = new TTree("SkimmingList", "SkimmingList");
	fSkimmingList->Branch("FileName", &fFileName);
	fSkimmingList->Branch("EventNumberInFile", &fEventNumberInFile);
	fSkimmingList->Branch("RunNumber", &fRunNumber);
	fSkimmingList->Branch("TimeStamp", &fEventTime);
	PostData(3, fSkimmingList);

	// Get input information
	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	fInputTree = am->GetInputEventHandler()->GetTree();
	if (!fInputTree) {
		puts("AliAnalysisTaskCDskimESD: input tree not found\n");
		return;
	}
	fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");


	// prevent the task from being run on proof
	if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
	    AliAnalysisManager::kProofAnalysis) {
		AliFatal("AliAnalysisTaskCDskimESD: cannot be run on PROOF!");
	}

	TFile *file = 0x0;
	file = OpenFile(1); // open file for the first output slot
	if (!file) {
		puts("AliAnalysisTaskCDskimESD: file not opened\n");
			return;
	}

	fOutputTree = fInputTreeCopy->CloneTree(0);
	if (!fOutputTree) {
		puts("AliAnalysisTaskCDskimESD: cloning tree not successful\n");
		return;
	}

	file->SetCompressionLevel(5); // caused a crash
	fOutputTree->SetDirectory(file);
	fOutputTree->SetAutoFlush(-10*1024*1024);

	if (fDoGapCond) {
		if(!fOutputTree->GetBranch("gapCondition")) {
			fOutputTree->Branch("gapCondition", &fCurrentGapCondition);
		}
	}

	PostData(1, fOutputTree);

	// process the mc information
	if(fDoMC) {
		fWrkDir = gDirectory; // save the old directory
		fKinematicsFile = new TFile("Kinematics.root", "update");
		gDirectory = fWrkDir; // restore the old directory

		fTrackRefsFile = new TFile("TrackRefs.root", "update");
		gDirectory = fWrkDir; // restore the old directory

		fTEinput =
			((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
		if (!fTEinput) {
			// this trick is done in order to post an output tree, although the tree
			// which is cloned is not available during the UserCreateOutputObjects()
			fTE = new TTree("TE", "TE");
		}
		else {
			fTE = fTEinput->CloneTree(0);
			//CopygALICE(am); // not needed, not properly working/tested
		}
		fgaliceFile = OpenFile(4);
		if (!file) {
			puts("AliAnalysisTaskCDskimESD: galice.root file not opened\n");
			return;
		}
		fgaliceFile->SetCompressionLevel(5);
		fTE->SetDirectory(fgaliceFile);
		fTE->SetAutoFlush(-10*1024*1024);


		PostData(4, fTE);
	}
}


//______________________________________________________________________________
/* // not needed => not yet working / properly tested
void AliAnalysisTaskCDskimESD::CopygALICE(AliAnalysisManager* am)
{
	//
	// copy all objects contained in the galice root excluding the TE tree
	//

	TString galiceInput =
		*(((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetInputPath());
	galiceInput += "galice.root";
	printf("galiceInput=%s\n", galiceInput.Data());
	fWrkDir = gDirectory;
	TFile* input = new TFile(galiceInput.Data(), "READ");

	input->cd();
	gDirectory->ls();

	//loop on all entries of of the input directory
	TKey *key;
	TIter nextkey(input->GetListOfKeys());
	input->GetListOfKeys()->ls();
	while ((key = (TKey*)nextkey())) {
		TString name = key->GetName();
		printf("name=%s\n", name.Data());
		const TString classname = key->GetClassName();
		printf("classname=%s\n", classname.Data());
		TClass *cl = gROOT->GetClass(classname);
		if (!cl) continue;
		if (!name.Contains("TE") && !classname.Contains("AliRun")) {
			TObject *obj = key->ReadObj();
			fgaliceFile->cd();
			obj->Write();
			delete obj;
			input->cd();
		}
		else if (cl->InheritsFrom("AliRunLoader")) {
			// TODO write the runloader stuff
			//AliRunLoader* rl = AliRunLoader::Instance();
			//printf("FILENAME=%s\n", rl->GetFileName().Data());
		}
		fgaliceFile->ls();
	}
	fgaliceFile->SaveSelf(kTRUE);
	fWrkDir->cd();
}
*/
Commit	Line	Data
369562e9	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	// $Id$
	17	//
	18	// Task to skim ESD files, basic idea found at PWGGA/EMCALTasks as
	19	// AliEsdSkimTask. This version is modified to produce a exact copy of the ESD
	20	// stream containing only events which are of interest for the central
	21	// diffractive analysis (PWGUD) and a control sample. In order to be able to
	22	// process MC data as well, the usual concept using an output slot is for
	23	// output data is violated and an output file is directly created without
	24	// a corresponding output slot.
	25	//
	26	// Author:
	27	// Felix Reidt <Felix.Reidt@cern.ch>
	28
	29	// ROOT headers
	30	#include <TClonesArray.h>
	31	#include <TFile.h>
	32	#include <TTree.h>
	33	#include <TH1.h>
	34	#include <TAxis.h>
	35	#include <TRandom3.h>
	36	//#include <TKey.h>
	37	//#include <TROOT.h>
	38
	39	// AliRoot headers
	40	#include "AliAnalysisManager.h"
	41	#include "AliCentrality.h"
	42	#include "AliESDEvent.h"
	43	#include "AliESDFMD.h"
	44	#include "AliEventplane.h"
	45	#include "AliInputEventHandler.h"
	46	#include "AliMultiplicity.h"
	47	#include "AliVEventHandler.h"
	48	#include "AliMCEventHandler.h"
	49	//#include "AliMCEvent.h"
	50	//#include "AliStack.h"
	51	//#include "AliGenEventHeader.h"
	52	//#include "AliRunLoader.h"
	53
	54	// own classes
	55	#include "AliCDMesonBase.h"
	56	#include "AliCDMesonUtils.h"
	57
	58	// own header
	59	#include "AliAnalysisTaskCDskimESD.h"
	60
	61	ClassImp(AliAnalysisTaskCDskimESD)
	62
	63	//______________________________________________________________________________
	64	AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD(const char *opt,
65	Bool_t reduceGapEvents)
66	: AliAnalysisTaskSE(opt)
67	, fInputTree(0x0)
68	, fInputTreeCopy(0x0)
69	, fOutputTree(0x0)
70	, fDoMC(kFALSE)
71	, fWrkDir(0x0)
72	, fKinematicsFile(0x0)
73	, fTrackRefsFile(0x0)
74	, fTEinput(0x0)
75	, fTE(0x0)
76	, fgaliceFile(0x0)
77	, fDoGapCond(kTRUE)
78	, fCurrentGapCondition(0)
79	, fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
80	// only events containing a V0 double gap or
81	, fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
82	, fReduceGapEvents(reduceGapEvents)
83	, fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
84	AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
85	, fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
86	, fStatsFlow(0x0)
87	, fSkimmingList(0x0)
88	, fFileName()
89	, fEventNumberInFile(-1)
90	, fRunNumber(-1)
91	, fEventTime(0)
92	{
93	//
94	// Constructor.
95	//
96
97	if (!opt)
98	return;
99
100	//fBranchNames = "ESD:AliESDHeader.,AliESDRun."; // TODO don't we need all?
101
102	// check whether we are running on MC data
103	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
104	AliMCEventHandler* mcH =
105	dynamic_cast<AliMCEventHandler*>(am->GetMCtruthEventHandler());
106	if(mcH) {
107	fDoMC = kTRUE;
108
109	fWrkDir = gDirectory; // save the old directory
110	fKinematicsFile = new TFile("Kinematics.root", "recreate");
111	fKinematicsFile->Close();
112	delete fKinematicsFile;
113	fKinematicsFile = 0x0;
114	gDirectory = fWrkDir; // restore the old directory
115
116	fTrackRefsFile = new TFile("TrackRefs.root", "recreate");
117	fTrackRefsFile->Close();
118	delete fTrackRefsFile;
119	fTrackRefsFile = 0x0;
120	gDirectory = fWrkDir; // restore the old directory
121	}
122
123	DefineOutput(1, TTree::Class());
124	DefineOutput(2, TH1::Class());
125	DefineOutput(3, TTree::Class());
126	if (fDoMC) DefineOutput(4, TTree::Class());
127	}
128
129
130	//______________________________________________________________________________
131	AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD()
132	: AliAnalysisTaskSE()
133	, fInputTree(0x0)
134	, fInputTreeCopy(0x0)
135	, fOutputTree(0x0)
136	, fDoMC(kFALSE)
137	, fWrkDir(0x0)
138	, fKinematicsFile(0x0)
139	, fTrackRefsFile(0x0)
140	, fTEinput(0x0)
141	, fTE(0x0)
142	, fgaliceFile(0x0)
143	, fDoGapCond(kTRUE)
144	, fCurrentGapCondition(0)
145	, fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
146	// only events containing a V0 double gap or
147	, fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
148	, fReduceGapEvents(kFALSE)
149	, fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
150	AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
151	, fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
152	, fStatsFlow(0x0)
153	, fSkimmingList(0x0)
154	, fFileName()
155	, fEventNumberInFile(-1)
156	, fRunNumber(-1)
157	, fEventTime(0)
158	{
159	//
160	// Default Constructor
161	//
162	}
163
164
165	//______________________________________________________________________________
166	AliAnalysisTaskCDskimESD::~AliAnalysisTaskCDskimESD()
167	{
168	//
169	// deconstructor
170	//
171
172	/* // this delete lead to chrases caused by the garbage collection ...
173	if (fOutputTree) {
174	delete fOutputTree;
175	fOutputTree = 0x0;
176	}
177	if (fInputTreeCopy) {
178	delete fInputTreeCopy;
179	fInputTreeCopy = 0x0;
180	}
181	if (fStatsFlow) {
182	delete fStatsFlow;
183	fStatsFlow = 0x0;
184	} */
185	if (fDoMC && fKinematicsFile) {
186	fKinematicsFile->Close();
187	delete fKinematicsFile;
188	fKinematicsFile = 0x0;
189	}
190	if (fDoMC && fTrackRefsFile) {
191	fTrackRefsFile->Close();
192	delete fTrackRefsFile;
193	fTrackRefsFile = 0x0;
194	}
195	}
196
197
198	//______________________________________________________________________________
199	void AliAnalysisTaskCDskimESD::UserExec(Option_t /opt*/)
200	{
201	//
202	// Process event.
203	//
204
205	fStatsFlow->Fill(0.); // UserExec(...) <= bin name
206
207	if (fDoMC) {
208	// replace dummy tree by the correct tree
209	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
210	if(!fTEinput) {
211	fTEinput =
212	((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
213	if (fTEinput) {
214	delete fTE;
215	fTE = fTEinput->CloneTree(0);
216	fTE->SetDirectory(fgaliceFile);
217	fTE->SetAutoFlush(-1010241024);
218	//CopygALICE(am);
219	}
220	}
221	else if (fTEinput != am->GetMCtruthEventHandler()->GetTree()) {
222	fTEinput =
223	((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
224	fTEinput->CopyAddresses(fTE);
225	}
226	}
227
228	// check whether the ESD input tree is still the same - NOT NECESSARY
229	//TTree* currentInput =
230	// AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree();
231	//if (fInputTree != currentInput) { // ESD input tree has changed
232	// puts("CHANGED INPUT TREE");
233	// fInputTree = currentInput;
234	// fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");
235	// fInputTreeCopy->CopyAddresses(fOutputTree);
236	//}
237
238	AliESDEvent esdin = dynamic_cast<AliESDEvent>(InputEvent());
b24ae32a	239	if (!esdin \|\| !fInputTree \|\| !fInputTreeCopy \|\| !fOutputTree \|\|
b24ae32a	240	(fDoMC && (!fTEinput \|\| !fTE))) {
369562e9	241	// check whether everything is ready
	242	PostData(1, fOutputTree);
	243	PostData(2, fStatsFlow);
	244	PostData(3, fSkimmingList);
	245	if (fDoMC) PostData(4, fTE);
	246	return;
	247	}
	248
	249	fStatsFlow->Fill(1.); // Ready <= bin name
	250
	251	//============================================================================
	252	// event selection, every event passing the next lines is written to the new
	253	// ESD stream
	254	//
	255	// modify this lines in order to adjust the ESD selection to your needs
	256
	257	// determine current gap condition
	258	fCurrentGapCondition = AliCDMesonUtils::GetGapConfig(esdin, 0x0, 0x0, 0x0,
	259	0x0, 0x0, 0x0, 0x0, 0x0);
	260
	261	TRandom3 rnd(0);
	262	if (!(AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) \|\|
	263	(fRequestedGapCondition & fCurrentGapCondition) \|\| fReduceGapEvents) {
	264	// check whether the event satifies the required gap condition and whether
	265	// all double-gap events will be stored
	266
	267	if (fReduceGapEvents &&
	268	(AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) &&
	269	!(fRequestedGapCondition & fCurrentGapCondition)) {
	270	// not all double gap events can be stored, but the gap determination
	271	// was successful and at least a according to the less stringent condition
	272	// a double gap was found
	273	if (!(fRefinedGapCondition & fCurrentGapCondition)) {
	274	// event fulfilled the more stringent gap and hence is will be stored
	275	fStatsFlow->Fill(5.);
	276	}
	277	else if (rnd.Rndm() > (1.-fReducedGapFraction)) {
	278	// randomly select a fraction of
	279	fStatsFlow->Fill(4.); // Double Gap Event <= bin name
	280	}
	281	else {
	282	// reject event (no refined gap and not randomly selected)
	283	fStatsFlow->Fill(3.); // Event Rejected <= bin name
	284	PostData(1, fOutputTree);
	285	PostData(2, fStatsFlow);
	286	PostData(3, fSkimmingList);
	287	if (fDoMC) PostData(4, fTE);
	288	return;
	289	}
	290	}
	291	else if (rnd.Rndm() > (1.-fNoGapFraction)) { // randomly selected
	292	fStatsFlow->Fill(2.); // Control Sample Event <= bin name
	293	}
	294	else {
	295	fStatsFlow->Fill(3.); // Event Rejected <= bin name
	296	PostData(1, fOutputTree);
	297	PostData(2, fStatsFlow);
	298	PostData(3, fSkimmingList);
	299	if (fDoMC) PostData(4, fTE);
	300	return;
	301	}
	302	}
	303	else {
	304	fStatsFlow->Fill(4.); // Double Gap Event <= bin name
305	}
306	// end of event selection
307	//============================================================================
308
309	// load the current event in the cloned input tree
310	//
311	// unfortunately the Entry() gives the entry number with in the current input
312	// tree, not within the chain for MC data
313	// => use GetReadEntry() of the ESD tree (value return by Entry() for real
314	//data)
315	fInputTreeCopy->GetEvent(fInputTree->GetReadEntry());
316	fOutputTree->Fill(); // fill the current event into the new ESD stream
317	//printf("Entry()=%lld\nEntries=%lld\nChainOffset=%lld\nReadEntry=%lld\n",
318	// Entry(), fInputTreeCopy->GetEntries(),
319	// fInputTreeCopy->GetChainOffset(), fInputTree->GetReadEntry());
320
321	// MC specific stuff ---------------------------------------------------------
322	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
323
324	if (fDoMC) { // do the MC related part of the skimming
325	Long64_t iEntry = fOutputTree->GetEntries() - 1;
326	AliMCEventHandler* mcHandler =
327	(AliMCEventHandler*)am->GetMCtruthEventHandler();
328	TTree* treeK = mcHandler->TreeK();
329	TTree* treeTR = mcHandler->TreeTR();
7d348e8a	330	if (!treeK) {
	331	AliFatal("TreeK not found!");
	332	return;
	333	}
	334	if (!treeTR) {
	335	AliFatal("TreeTR not found!");
	336	return;
	337	}
369562e9	338
	339	fWrkDir = gDirectory;
	340	gDirectory = fKinematicsFile;
	341	TString dirName = TString(Form("Event%ld", (Long_t)iEntry));
	342	gDirectory->mkdir(dirName);
	343	gDirectory->cd(dirName);
	344	TTree* outputTreeK = treeK->CloneTree(); // copy the whole tree
7d348e8a	345	if (!outputTreeK) {
	346	AliFatal("Could not clone TreeK!");
	347	return;
	348	}
369562e9	349	outputTreeK->Write();
	350	treeK->CopyAddresses(outputTreeK, kTRUE); // separate clone again
	351	treeK->GetListOfClones()->Remove((TObject*)outputTreeK);
	352	outputTreeK->Delete();
	353	outputTreeK = 0x0;
	354
	355	gDirectory = fTrackRefsFile;
	356	gDirectory->mkdir(dirName);
	357	gDirectory->cd(dirName);
	358	TTree* outputTreeTR = treeTR->CloneTree(); // copy the whole tree
7d348e8a	359	if (!outputTreeTR) {
	360	AliFatal("Could not clone TreeTR");
	361	return;
	362	}
369562e9	363	outputTreeTR->Write();
	364	treeTR->CopyAddresses(outputTreeTR, kTRUE); // separate clone again
	365	treeTR->GetListOfClones()->Remove((TObject*)outputTreeTR);
	366	outputTreeTR->Delete();
	367	outputTreeTR = 0x0;
	368	gDirectory = fWrkDir;
	369
	370	TTree* currentTEinput =
	371	((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
	372	if (fTEinput != currentTEinput) {
	373	// when the input file is changed, the kinematics tree address changes
	374	// hence our output tree fTE has to be linked to the new one!
	375	fTEinput = currentTEinput;
	376	fTEinput->CopyAddresses(fTE);
	377	}
	378	fTE->Fill(); // fill MC event headers
	379	}
	380
	381	fStatsFlow->Fill(6.); // Tree Filled <= bin name
	382
	383	// remember which event we stored
	384	fEventNumberInFile = esdin->GetEventNumberInFile();
	385	fRunNumber = esdin->GetRunNumber();
	386	fFileName = fInputTree->GetCurrentFile()->GetName();
	387	fEventTime = esdin->GetTimeStamp();
	388	fSkimmingList->Fill();
	389
	390	PostData(1, fOutputTree);
	391	PostData(2, fStatsFlow);
	392	PostData(3, fSkimmingList);
	393	if (fDoMC) PostData(4, fTE);
	394	}
	395
	396	//______________________________________________________________________________
	397	void AliAnalysisTaskCDskimESD::UserCreateOutputObjects()
	398	{
	399	// Create output objects.
	400	fStatsFlow = (TH1*)(new TH1F("skimmingStatsFlow", "", 7, 0., 7.));
	401	TAxis* x = fStatsFlow->GetXaxis();
	402	x->SetBinLabel(1, "UserExec(...)");
	403	x->SetBinLabel(2, "Ready");
	404	x->SetBinLabel(3, "Control Sample Event");
	405	x->SetBinLabel(4, "Rejected Event");
	406	x->SetBinLabel(5, "Double Gap Event");
	407	x->SetBinLabel(6, "Refined Double Gap Event");
	408	x->SetBinLabel(7, "Tree Filled");
	409
	410	PostData(2, fStatsFlow);
	411
	412	// TODO implement fSkimmingList
	413	fSkimmingList = new TTree("SkimmingList", "SkimmingList");
	414	fSkimmingList->Branch("FileName", &fFileName);
	415	fSkimmingList->Branch("EventNumberInFile", &fEventNumberInFile);
	416	fSkimmingList->Branch("RunNumber", &fRunNumber);
	417	fSkimmingList->Branch("TimeStamp", &fEventTime);
	418	PostData(3, fSkimmingList);
	419
	420	// Get input information
	421	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	422	fInputTree = am->GetInputEventHandler()->GetTree();
	423	if (!fInputTree) {
	424	puts("AliAnalysisTaskCDskimESD: input tree not found\n");
	425	return;
	426	}
427	fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");
428
429
430	// prevent the task from being run on proof
431	if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
432	AliAnalysisManager::kProofAnalysis) {
433	AliFatal("AliAnalysisTaskCDskimESD: cannot be run on PROOF!");
434	}
435
436	TFile *file = 0x0;
437	file = OpenFile(1); // open file for the first output slot
438	if (!file) {
439	puts("AliAnalysisTaskCDskimESD: file not opened\n");
440	return;
441	}
442
443	fOutputTree = fInputTreeCopy->CloneTree(0);
444	if (!fOutputTree) {
445	puts("AliAnalysisTaskCDskimESD: cloning tree not successful\n");
446	return;
447	}
448
449	file->SetCompressionLevel(5); // caused a crash
450	fOutputTree->SetDirectory(file);
451	fOutputTree->SetAutoFlush(-1010241024);
452
453	if (fDoGapCond) {
454	if(!fOutputTree->GetBranch("gapCondition")) {
455	fOutputTree->Branch("gapCondition", &fCurrentGapCondition);
456	}
457	}
458
459	PostData(1, fOutputTree);
460
461	// process the mc information
462	if(fDoMC) {
463	fWrkDir = gDirectory; // save the old directory
464	fKinematicsFile = new TFile("Kinematics.root", "update");
465	gDirectory = fWrkDir; // restore the old directory
466
467	fTrackRefsFile = new TFile("TrackRefs.root", "update");
468	gDirectory = fWrkDir; // restore the old directory
469
470	fTEinput =
471	((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
472	if (!fTEinput) {
473	// this trick is done in order to post an output tree, although the tree
474	// which is cloned is not available during the UserCreateOutputObjects()
475	fTE = new TTree("TE", "TE");
476	}
477	else {
478	fTE = fTEinput->CloneTree(0);
479	//CopygALICE(am); // not needed, not properly working/tested
480	}
481	fgaliceFile = OpenFile(4);
482	if (!file) {
483	puts("AliAnalysisTaskCDskimESD: galice.root file not opened\n");
484	return;
485	}
486	fgaliceFile->SetCompressionLevel(5);
487	fTE->SetDirectory(fgaliceFile);
488	fTE->SetAutoFlush(-1010241024);
489
490
491	PostData(4, fTE);
492	}
493	}
494
495
496	//______________________________________________________________________________
497	/* // not needed => not yet working / properly tested
498	void AliAnalysisTaskCDskimESD::CopygALICE(AliAnalysisManager* am)
499	{
500	//
501	// copy all objects contained in the galice root excluding the TE tree
502	//
503
504	TString galiceInput =
505	(((AliMCEventHandler)am->GetMCtruthEventHandler())->GetInputPath());
506	galiceInput += "galice.root";
507	printf("galiceInput=%s\n", galiceInput.Data());
508	fWrkDir = gDirectory;
509	TFile* input = new TFile(galiceInput.Data(), "READ");
510
511	input->cd();
512	gDirectory->ls();
513
514	//loop on all entries of of the input directory
515	TKey *key;
516	TIter nextkey(input->GetListOfKeys());
517	input->GetListOfKeys()->ls();
518	while ((key = (TKey*)nextkey())) {
519	TString name = key->GetName();
520	printf("name=%s\n", name.Data());
521	const TString classname = key->GetClassName();
522	printf("classname=%s\n", classname.Data());
523	TClass *cl = gROOT->GetClass(classname);
524	if (!cl) continue;
525	if (!name.Contains("TE") && !classname.Contains("AliRun")) {
526	TObject *obj = key->ReadObj();
527	fgaliceFile->cd();
528	obj->Write();
529	delete obj;
530	input->cd();
531	}
532	else if (cl->InheritsFrom("AliRunLoader")) {
533	// TODO write the runloader stuff
534	//AliRunLoader* rl = AliRunLoader::Instance();
535	//printf("FILENAME=%s\n", rl->GetFileName().Data());
536	}
537	fgaliceFile->ls();
538	}
539	fgaliceFile->SaveSelf(kTRUE);
540	fWrkDir->cd();
541	}
542	*/