[u/mrichter/AliRoot.git] / PWG0 / multVScentPbPb / MyAnalysisMacro.C

void MyAnalysisMacro(TString dataset="/alice/sim/LHC10f8f_130844",
		     TString outFName="trbg.root",
		     Bool_t doRec  = kTRUE,
		     Bool_t doInj  = kTRUE,
		     Bool_t doRot  = kTRUE,
		     Bool_t doMix  = kFALSE,
		     Bool_t useMC  = kTRUE,
		     Bool_t checkReconstructables = kTRUE,
		     // 
		     Float_t etaCut     = 2.5,
		     //
		     // specific parameters for reconstruction
		     //----------------------- Zv selection parameters important for mixing, to be tuned
		     Float_t zMin       = -20,
		     Float_t zMax       =  20,
		     Float_t zMixBinSz  =  5, //0.1,
		     //---------------------------------------------------------------------------------
		     //
		     //----------------------- Ntracklets selection parameters important for mixing, to be tuned
		     Float_t ntMin      =   1,
		     Float_t ntMax      = 15000,
		     Float_t ntMixBinSz =  5000, 
		     //---------------------------------------------------------------------------------
		     //
		     float  phiRot      = 3.14159e+00,
		     float  injScale    = 0.7,
		     Bool_t scaleDTheta = kTRUE,
		     float  nStdDev     = 25.,
		     float  dphi        = 0.08, 
		     float  dtht        = 0.025, 
		     float  phishift    = 0.0045, 
		     Bool_t remOvl      = kTRUE, 
		     float  ovlPhiCut   = 0.005, 
		     float  ovlZetaCut  = 0.05,
		     Int_t  nEvents     = 50000,
		     Int_t  nEventsSkip = 0) 
{
  //  
  TString format = GetFormatFromDataSet(dataset);
  //
  // ALICE stuff
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  if (!mgr) mgr = new AliAnalysisManager("Test train");
  //
  InputHandlerSetup(format,useMC);
  if (doMix) MixHandlerSetup(ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
  // compile our task
  gProof->Load("AliITSMultRecBg.cxx++");
  gProof->Load("AliTrackletTaskUni.cxx++");
  //
  // load and run AddTask macro
  gROOT->LoadMacro("AddMultTaskRS.C");
  //
  // create our task
  AliTrackletTaskUni *mltTask = AddMultTaskRS(outFName.Data());
  //
  mltTask->SetDoNormalReco(doRec);
  mltTask->SetDoInjection(doInj);
  mltTask->SetDoRotation(doRot);
  mltTask->SetDoMixing(doMix);  
  //
  mltTask->SetUseMC(useMC);
  mltTask->SetCheckReconstructables(checkReconstructables);
  //
  mltTask->SetEtaCut(etaCut);
  mltTask->SetZVertexMin(zMin);
  mltTask->SetZVertexMax(zMax);
  mltTask->SetMultCutMin(ntMin);
  mltTask->SetMultCutMax(ntMax);
  //
  mltTask->SetScaleDThetaBySin2T(scaleDTheta);
  mltTask->SetNStdDev(nStdDev);
  mltTask->SetPhiWindow(dphi);
  mltTask->SetThetaWindow(dtht);
  mltTask->SetPhiShift(phishift);
  mltTask->SetPhiOverlapCut(ovlPhiCut);
  mltTask->SetZetaOverlapCut(ovlZetaCut);
  mltTask->SetPhiRot(phiRot);
  mltTask->SetInjScale(injScale);
  mltTask->SetRemoveOverlaps(remOvl);
  //
  printf("new Task: %p\n",mltTask);
  //
  AddPhysicsSelection(useMC);
  mltTask->SelectCollisionCandidates(useMC ? (AliVEvent::kMB) : (AliVEvent::kUserDefined) );
  //
  // Run analysis
  mgr->InitAnalysis();
  // process dataset  
  mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip); 
  //
  TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_"; 
  evstCmd += outFName;  evstCmd += " \nfi";
  gSystem->Exec( evstCmd.Data() );
  
}


TString GetFormatFromDataSet(TString dataset) {
  
//   Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
  TString dsTreeName;
  if (dataset.Contains("#")) {
    Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
    dsTreeName=dataset(dataset.Last('#'),dataset.Length());
  } else {
    Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
    TFileCollection *ds = gProof->GetDataSet(dataset.Data());
    if (!ds) {
      Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
      return "";
    }
    dsTreeName = ds->GetDefaultTreeName();
  }

  if (dsTreeName.Contains("esdTree")) {
    Info("runAAF.C","ESD input format detected ...");
    return "ESD";
  } else if (dsTreeName.Contains("aodTree"))  {
    Info("runAAF.C","AOD input format detected ...");
    return "AOD";
  } else {
    Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
    Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
  }
  
  return "";
}

Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
{
  format.ToLower();

  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();

  AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();

  if (cin) return;

  if (!format.CompareTo("esd"))
  {
    AliESDInputHandler *esdInputHandler = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

    if (!esdInputHandler)
    {
      Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
      esdInputHandler = new AliESDInputHandlerRP();
      mgr->SetInputEventHandler(esdInputHandler);
    }

    if (useKine)
    {
      AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());

      if (!mcInputHandler)
      {
        Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
        AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
        mgr->SetMCtruthEventHandler(mcInputHandler);
      }
    }

  }
  else if (!format.CompareTo("aod"))
  {
    AliAODInputHandler *aodInputHandler = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

    if (!aodInputHandler)
    {
      Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
      aodInputHandler = new AliAODInputHandler();
      mgr->SetInputEventHandler(aodInputHandler);
    }
  }
  else
  {
    Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
    return kFALSE;
  }

  return kTRUE;
}

void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
		     float zMin, float zMax, float zMixBinSz)
{
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  if (!mgr) return;
  int bufferSize = 1;
  AliESDInputHandlerRP *esdH = dynamic_cast<AliESDInputHandlerRP*>(mgr->GetInputEventHandler());
  if (!esdH) return;
  //
  AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
  esdMixH->SetInputHandlerForMixing(esdH);
  AliMixEventPool *evPool = new AliMixEventPool("MyPool");
  AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
  AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
  //  evPool->AddCut(tracklets);
  evPool->AddCut(zvertex);
  //evPool->Init();
  evPool->Print();
  esdMixH->SetEventPool(evPool);
  esdH->SetMixingHandler(esdMixH);
}

void AddPhysicsSelection(Bool_t isMC)
{
  // physics selection a la Michele
  printf("Requesting physics selection in %s mode\n",isMC ? "MC":"Data");
  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
  //isMC is true when processing monte carlo, the second 0 disables the cluster vs tracklets
  AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(isMC,0);
  if(!isMC) {
    AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
    physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
    /*
    physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
    physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
    */
    //
    physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
    physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
    //
    // This are needed only to fill the statistics tables
    physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
    //
    /*
    physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
    //
    physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
    //
    */
    /*
    //
    physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
    //
    physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
    physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
    */
  } 
  // if you use the following line, your task only gets the selected events
  //  task->SelectCollisionCandidates(AliVEvent::kUserDefined);
  //
  //Alternatively, in the UserExec of your task:
  //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
  //
}
Commit	Line	Data
a9a39f46	1	void MyAnalysisMacro(TString dataset="/alice/sim/LHC10f8f_130844",
	2	TString outFName="trbg.root",
	3	Bool_t doRec = kTRUE,
	4	Bool_t doInj = kTRUE,
	5	Bool_t doRot = kTRUE,
	6	Bool_t doMix = kFALSE,
	7	Bool_t useMC = kTRUE,
	8	Bool_t checkReconstructables = kTRUE,
	9	//
	10	Float_t etaCut = 2.5,
	11	//
	12	// specific parameters for reconstruction
	13	//----------------------- Zv selection parameters important for mixing, to be tuned
	14	Float_t zMin = -20,
	15	Float_t zMax = 20,
	16	Float_t zMixBinSz = 5, //0.1,
	17	//---------------------------------------------------------------------------------
	18	//
	19	//----------------------- Ntracklets selection parameters important for mixing, to be tuned
	20	Float_t ntMin = 1,
	21	Float_t ntMax = 15000,
	22	Float_t ntMixBinSz = 5000,
	23	//---------------------------------------------------------------------------------
	24	//
	25	float phiRot = 3.14159e+00,
	26	float injScale = 0.7,
	27	Bool_t scaleDTheta = kTRUE,
	28	float nStdDev = 25.,
	29	float dphi = 0.08,
	30	float dtht = 0.025,
	31	float phishift = 0.0045,
	32	Bool_t remOvl = kTRUE,
	33	float ovlPhiCut = 0.005,
	34	float ovlZetaCut = 0.05,
	35	Int_t nEvents = 50000,
	36	Int_t nEventsSkip = 0)
	37	{
	38	//
	39	TString format = GetFormatFromDataSet(dataset);
	40	//
	41	// ALICE stuff
	42	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	43	if (!mgr) mgr = new AliAnalysisManager("Test train");
	44	//
	45	InputHandlerSetup(format,useMC);
	46	if (doMix) MixHandlerSetup(ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
	47	// compile our task
	48	gProof->Load("AliITSMultRecBg.cxx++");
	49	gProof->Load("AliTrackletTaskUni.cxx++");
	50	//
	51	// load and run AddTask macro
	52	gROOT->LoadMacro("AddMultTaskRS.C");
	53	//
	54	// create our task
	55	AliTrackletTaskUni *mltTask = AddMultTaskRS(outFName.Data());
	56	//
	57	mltTask->SetDoNormalReco(doRec);
	58	mltTask->SetDoInjection(doInj);
	59	mltTask->SetDoRotation(doRot);
	60	mltTask->SetDoMixing(doMix);
	61	//
	62	mltTask->SetUseMC(useMC);
	63	mltTask->SetCheckReconstructables(checkReconstructables);
	64	//
65	mltTask->SetEtaCut(etaCut);
66	mltTask->SetZVertexMin(zMin);
67	mltTask->SetZVertexMax(zMax);
68	mltTask->SetMultCutMin(ntMin);
69	mltTask->SetMultCutMax(ntMax);
70	//
71	mltTask->SetScaleDThetaBySin2T(scaleDTheta);
72	mltTask->SetNStdDev(nStdDev);
73	mltTask->SetPhiWindow(dphi);
74	mltTask->SetThetaWindow(dtht);
75	mltTask->SetPhiShift(phishift);
76	mltTask->SetPhiOverlapCut(ovlPhiCut);
77	mltTask->SetZetaOverlapCut(ovlZetaCut);
78	mltTask->SetPhiRot(phiRot);
79	mltTask->SetInjScale(injScale);
80	mltTask->SetRemoveOverlaps(remOvl);
81	//
82	printf("new Task: %p\n",mltTask);
83	//
84	AddPhysicsSelection(useMC);
85	mltTask->SelectCollisionCandidates(useMC ? (AliVEvent::kMB) : (AliVEvent::kUserDefined) );
86	//
87	// Run analysis
88	mgr->InitAnalysis();
89	// process dataset
90	mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip);
91	//
92	TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_";
93	evstCmd += outFName; evstCmd += " \nfi";
94	gSystem->Exec( evstCmd.Data() );
95
96	}
97
98
99	TString GetFormatFromDataSet(TString dataset) {
100
101	// Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
102	TString dsTreeName;
103	if (dataset.Contains("#")) {
104	Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
105	dsTreeName=dataset(dataset.Last('#'),dataset.Length());
106	} else {
107	Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
108	TFileCollection *ds = gProof->GetDataSet(dataset.Data());
109	if (!ds) {
110	Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
111	return "";
112	}
113	dsTreeName = ds->GetDefaultTreeName();
114	}
115
116	if (dsTreeName.Contains("esdTree")) {
117	Info("runAAF.C","ESD input format detected ...");
118	return "ESD";
119	} else if (dsTreeName.Contains("aodTree")) {
120	Info("runAAF.C","AOD input format detected ...");
121	return "AOD";
122	} else {
123	Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
124	Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
125	}
126
127	return "";
128	}
129
130	Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
131	{
132	format.ToLower();
133
134	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
135
136	AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
137
138	if (cin) return;
139
140	if (!format.CompareTo("esd"))
141	{
142	AliESDInputHandler esdInputHandler = dynamic_cast<AliESDInputHandler>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
143
144	if (!esdInputHandler)
145	{
146	Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
147	esdInputHandler = new AliESDInputHandlerRP();
148	mgr->SetInputEventHandler(esdInputHandler);
149	}
150
151	if (useKine)
152	{
153	AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
154
155	if (!mcInputHandler)
156	{
157	Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
158	AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
159	mgr->SetMCtruthEventHandler(mcInputHandler);
160	}
161	}
162
163	}
164	else if (!format.CompareTo("aod"))
165	{
166	AliAODInputHandler aodInputHandler = dynamic_cast<AliAODInputHandler>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
167
168	if (!aodInputHandler)
169	{
170	Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
171	aodInputHandler = new AliAODInputHandler();
172	mgr->SetInputEventHandler(aodInputHandler);
173	}
174	}
175	else
176	{
177	Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
178	return kFALSE;
179	}
180
181	return kTRUE;
182	}
183
184	void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
185	float zMin, float zMax, float zMixBinSz)
186	{
187	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
188	if (!mgr) return;
189	int bufferSize = 1;
190	AliESDInputHandlerRP esdH = dynamic_cast<AliESDInputHandlerRP>(mgr->GetInputEventHandler());
191	if (!esdH) return;
192	//
193	AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
194	esdMixH->SetInputHandlerForMixing(esdH);
195	AliMixEventPool *evPool = new AliMixEventPool("MyPool");
196	AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
197	AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
198	// evPool->AddCut(tracklets);
199	evPool->AddCut(zvertex);
200	//evPool->Init();
201	evPool->Print();
202	esdMixH->SetEventPool(evPool);
203	esdH->SetMixingHandler(esdMixH);
204	}
205
206	void AddPhysicsSelection(Bool_t isMC)
207	{
208	// physics selection a la Michele
209	printf("Requesting physics selection in %s mode\n",isMC ? "MC":"Data");
210	gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
211	//isMC is true when processing monte carlo, the second 0 disables the cluster vs tracklets
212	AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(isMC,0);
213	if(!isMC) {
214	AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
215	physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
216	/*
217	physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
218	physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
219	*/
220	//
221	physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
222	physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
223	//
224	// This are needed only to fill the statistics tables
225	physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
226	physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
227	physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
228	//
229	/*
230	physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
231	physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
232	physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
233	//
234	physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
235	physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
236	physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
237	//
238	*/
239	/*
240	//
241	physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
242	physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
243	physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
244	//
245	physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
246	physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
247	physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
248	*/
249	}
250	// if you use the following line, your task only gets the selected events
251	// task->SelectCollisionCandidates(AliVEvent::kUserDefined);
252	//
253	//Alternatively, in the UserExec of your task:
254	//Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
255	//
256	}