[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEFluctuationAnalysisTaskTrain.cxx

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TParticle.h"
#include "TObjArray.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
//#include "AliESDtrackCuts.h"
#include "AliCentrality.h"

#include "AliEbyEFluctuationAnalysisTaskTrain.h"

// Event by event charge fluctuation analysis
// Authors: Satyajit Jena and Panos Cristakoglou
// 

ClassImp(AliEbyEFluctuationAnalysisTaskTrain)

//________________________________________________________________________
AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name) 
  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), 
    fHistEventStats(0), fHistCentrality(0),
    fHistNMult(0), fHistNPlusNMinus(0),
    fHistNMultMC(0), fHistNPlusNMinusMC(0),
    fAnalysisType("ESD"), fAnalysisMode("TPC"),
    fCentralityEstimator("V0M"), 
    fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),
    fVxMax(3.), fVyMax(3.), fVzMax(10.),
    fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), 
    fMaxDCAxy(3.0), fMaxDCAz(3.0) {
  // Constructor
  
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {
  // Create histograms
  // Called once

  fOutputList = new TList();
  fOutputList->SetOwner();

  //Event stats.
  TString gCutName[4] = {"Total","Offline trigger",
                         "Vertex","Analyzed"};
  fHistEventStats = new TH1F("fHistEventStats",
                             "Event statistics;;N_{events}",
                             4,0.5,4.5);
  for(Int_t i = 1; i <= 4; i++)
    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
  fOutputList->Add(fHistEventStats);

  //ESD analysis
  if(fAnalysisType.CompareTo("ESD") == 0 ) {
    fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
			       20,0.5,20.5);
    fOutputList->Add(fHistCentrality);
    
    TString histName;
    histName = "fHistNMult";
    fHistNMult = new TH1F(histName.Data(), 
			  ";N_{mult.}",
			  800,0,4000);
    fOutputList->Add(fHistNMult);

    histName = "fHistNPlusNMinus";
    fHistNPlusNMinus = new TH2F(histName.Data(), 
				";N_{+};N_{-}",
				2000,0.5,2000.5,2000,0.5,2000.5);  
    fOutputList->Add(fHistNPlusNMinus);
  }//ESD analysis
  else if(fAnalysisType.CompareTo("MC") == 0) {
    TString histName = "fHistNMultMC";
    fHistNMultMC = new TH1F(histName.Data(), 
			    ";N_{mult.}",
			    800,0,8000);
    fOutputList->Add(fHistNMultMC);
    
    histName = "fHistNPlusNMinusMC";
    fHistNPlusNMinusMC = new TH2F(histName.Data(), 
				  ";N_{+};N_{-}",
				  3000,0.5,3000.5,3000,0.5,3000.5);  
    fOutputList->Add(fHistNPlusNMinusMC);
  }//MC analysis

  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {
  // Main loop
  // Called for each event
  
  Int_t nPlus = 0, nMinus = 0;

  // Post output data.
  //ESD analysis
  if(fAnalysisType.CompareTo("ESD") == 0) {
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    if (!fESD) {
      printf("ERROR: fESD not available\n");
      return;
    }
  
    fHistEventStats->Fill(1); //all events
    
    //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
    //if(isSelected) {
      fHistEventStats->Fill(2); //triggered + centrality
      

      //Centrality stuff
      AliCentrality *centrality = fESD->GetCentrality();
      Int_t nCentrality = 0;
      
      if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
					      fCentralityPercentileMax,
					      fCentralityEstimator.Data())) {
	if(fAnalysisMode.CompareTo("TPC") == 0 ) {
	  const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
	  if(vertex) {
	    if(vertex->GetNContributors() > 0) {
	      if(vertex->GetZRes() != 0) {
		fHistEventStats->Fill(3); //events with a proper vertex
		if(TMath::Abs(vertex->GetX()) < fVxMax) {
		  if(TMath::Abs(vertex->GetY()) < fVyMax) {
		    if(TMath::Abs(vertex->GetZ()) < fVzMax) {
		      fHistEventStats->Fill(4); //analyzed events
		      
		      Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
                             centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
                             nCentrality,fESD->GetNumberOfTracks());
		      

		      Int_t nAcceptedTracks = 0;
		      
		      nAcceptedTracks = fESD->GetNumberOfTracks();
		      AliESDtrack* track = 0;
		      
		      Float_t dcaXY = 0.0, dcaZ = 0.0;
		      
		      for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) 
			{
			  track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));
			  if (!track) {
			    printf("ERROR: Could not receive track %d\n", iTracks);
			    continue;
			  }
			  
			  AliESDtrack *tpcOnlyTrack = new AliESDtrack();
			  
			  if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
			    delete tpcOnlyTrack;
			    continue;
			  }
			  
			  tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);
			  			  
			  Int_t nClustersITS = track->GetITSclusters(0x0);
			  Int_t nClustersTPC = track->GetTPCclusters(0x0);
			  
			  Float_t chi2PerClusterITS = -1;
			  if (nClustersITS!=0)
			    chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
			  Float_t chi2PerClusterTPC = -1;
			  if (nClustersTPC!=0)
			    chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
			  
			  //	
			  
			  if(nClustersTPC < fMinNumberOfTPCClusters) continue;
			  if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;
			  if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;
			  if(TMath::Abs(dcaZ) > fMaxDCAz) continue;
			  
			  Short_t gCharge = tpcOnlyTrack->Charge();
			  if(gCharge > 0) nPlus += 1;
			  if(gCharge < 0) nMinus += 1;
			  delete tpcOnlyTrack;
			} // track loop
		      
		      if( nPlus == 0 && nMinus == 0) return;

		      //  cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;
		      fHistCentrality->Fill(nCentrality);
		      fHistNPlusNMinus->Fill(nPlus,nMinus);
		      fHistNMult->Fill(nPlus+nMinus);
		      
		      
		    }//Vz cut
		  }//Vy cut
		}//Vx cut
	      }//Vz resolution
	    }//number of contributors
	  }//valid vertex
	}//TPC analysis mode
      }//centrality
      //}//physics selection
  }//ESD analysis level
  
      //MC analysis
  if(fAnalysisType.CompareTo("MC") == 0 ) {
    AliMCEvent* mcEvent = MCEvent();
    if (!mcEvent) {
      Printf("ERROR: Could not retrieve MC event");
      return;
    }
    AliStack* stack = mcEvent->Stack();
    if (!stack) {
      Printf("ERROR: Could not retrieve MC stack");
      return;
    }
    
    fHistEventStats->Fill(1); 
    fHistEventStats->Fill(2); 
    fHistEventStats->Fill(3); 
    fHistEventStats->Fill(4); //analyzed events
    for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
      AliVParticle* track = mcEvent->GetTrack(iTracks);
      if (!track) {
	Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
	continue;
      }
      
      if(!stack->IsPhysicalPrimary(iTracks)) continue;
      if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
      if(TMath::Abs(track->Eta()) > 0.8) continue;
      
      Short_t gCharge = track->Charge();
      if(gCharge > 0) nPlus += 1;
      if(gCharge < 0) nMinus += 1;
    }//particle loop
    fHistNPlusNMinusMC->Fill(nPlus,nMinus);
    fHistNMultMC->Fill(nPlus+nMinus);
    
  }//MC analysis level


}      

//________________________________________________________________________
void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {
  // Draw result to the screen
  // Called once at the end of the query

  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    printf("ERROR: Output list not available\n");
    return;
  }
}
Commit	Line	Data
a65a7e70	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TCanvas.h"
	6	#include "TParticle.h"
	7	#include "TObjArray.h"
	8
	9	#include "AliAnalysisTask.h"
	10	#include "AliAnalysisManager.h"
	11
	12	#include "AliStack.h"
	13	#include "AliMCEvent.h"
	14	#include "AliESDEvent.h"
	15	#include "AliESDInputHandler.h"
	16	//#include "AliESDtrackCuts.h"
	17	#include "AliCentrality.h"
	18
	19	#include "AliEbyEFluctuationAnalysisTaskTrain.h"
	20
	21	// Event by event charge fluctuation analysis
	22	// Authors: Satyajit Jena and Panos Cristakoglou
	23	//
	24
	25	ClassImp(AliEbyEFluctuationAnalysisTaskTrain)
	26
	27	//________________________________________________________________________
	28	AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name)
	29	: AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
	30	fHistEventStats(0), fHistCentrality(0),
	31	fHistNMult(0), fHistNPlusNMinus(0),
	32	fHistNMultMC(0), fHistNPlusNMinusMC(0),
	33	fAnalysisType("ESD"), fAnalysisMode("TPC"),
	34	fCentralityEstimator("V0M"),
	35	fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),
	36	fVxMax(3.), fVyMax(3.), fVzMax(10.),
	37	fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0),
	38	fMaxDCAxy(3.0), fMaxDCAz(3.0) {
	39	// Constructor
	40
	41	// Define input and output slots here
	42	// Input slot #0 works with a TChain
	43	DefineInput(0, TChain::Class());
	44	// Output slot #0 id reserved by the base class for AOD
	45	// Output slot #1 writes into a TH1 container
	46	DefineOutput(1, TList::Class());
	47	}
	48
	49	//________________________________________________________________________
	50	void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {
	51	// Create histograms
	52	// Called once
	53
	54	fOutputList = new TList();
	55	fOutputList->SetOwner();
	56
	57	//Event stats.
	58	TString gCutName[4] = {"Total","Offline trigger",
	59	"Vertex","Analyzed"};
	60	fHistEventStats = new TH1F("fHistEventStats",
	61	"Event statistics;;N_{events}",
	62	4,0.5,4.5);
	63	for(Int_t i = 1; i <= 4; i++)
	64	fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
65	fOutputList->Add(fHistEventStats);
66
67	//ESD analysis
68	if(fAnalysisType.CompareTo("ESD") == 0 ) {
69	fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
70	20,0.5,20.5);
71	fOutputList->Add(fHistCentrality);
72
73	TString histName;
74	histName = "fHistNMult";
75	fHistNMult = new TH1F(histName.Data(),
76	";N_{mult.}",
77	800,0,4000);
78	fOutputList->Add(fHistNMult);
79
80	histName = "fHistNPlusNMinus";
81	fHistNPlusNMinus = new TH2F(histName.Data(),
82	";N_{+};N_{-}",
83	2000,0.5,2000.5,2000,0.5,2000.5);
84	fOutputList->Add(fHistNPlusNMinus);
85	}//ESD analysis
86	else if(fAnalysisType.CompareTo("MC") == 0) {
87	TString histName = "fHistNMultMC";
88	fHistNMultMC = new TH1F(histName.Data(),
89	";N_{mult.}",
90	800,0,8000);
91	fOutputList->Add(fHistNMultMC);
92
93	histName = "fHistNPlusNMinusMC";
94	fHistNPlusNMinusMC = new TH2F(histName.Data(),
95	";N_{+};N_{-}",
96	3000,0.5,3000.5,3000,0.5,3000.5);
97	fOutputList->Add(fHistNPlusNMinusMC);
98	}//MC analysis
99
100	PostData(1, fOutputList);
101	}
102
103	//________________________________________________________________________
104	void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {
105	// Main loop
106	// Called for each event
107
108	Int_t nPlus = 0, nMinus = 0;
109
110	// Post output data.
111	//ESD analysis
112	if(fAnalysisType.CompareTo("ESD") == 0) {
113	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
114	if (!fESD) {
115	printf("ERROR: fESD not available\n");
116	return;
117	}
118
119	fHistEventStats->Fill(1); //all events
120
121	//Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
122	//if(isSelected) {
123	fHistEventStats->Fill(2); //triggered + centrality
124
125
126	//Centrality stuff
127	AliCentrality *centrality = fESD->GetCentrality();
128	Int_t nCentrality = 0;
129
130	if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
131	fCentralityPercentileMax,
132	fCentralityEstimator.Data())) {
133	if(fAnalysisMode.CompareTo("TPC") == 0 ) {
134	const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
135	if(vertex) {
136	if(vertex->GetNContributors() > 0) {
137	if(vertex->GetZRes() != 0) {
138	fHistEventStats->Fill(3); //events with a proper vertex
e690d4d0	139	if(TMath::Abs(vertex->GetX()) < fVxMax) {
	140	if(TMath::Abs(vertex->GetY()) < fVyMax) {
	141	if(TMath::Abs(vertex->GetZ()) < fVzMax) {
a65a7e70	142	fHistEventStats->Fill(4); //analyzed events
	143
	144	Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
	145	centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
	146	nCentrality,fESD->GetNumberOfTracks());
	147
	148
	149	Int_t nAcceptedTracks = 0;
	150
	151	nAcceptedTracks = fESD->GetNumberOfTracks();
	152	AliESDtrack* track = 0;
	153
	154	Float_t dcaXY = 0.0, dcaZ = 0.0;
	155
	156	for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++)
	157	{
	158	track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));
	159	if (!track) {
	160	printf("ERROR: Could not receive track %d\n", iTracks);
	161	continue;
	162	}
	163
	164	AliESDtrack *tpcOnlyTrack = new AliESDtrack();
	165
	166	if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
	167	delete tpcOnlyTrack;
	168	continue;
	169	}
	170
	171	tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);
	172
	173	Int_t nClustersITS = track->GetITSclusters(0x0);
	174	Int_t nClustersTPC = track->GetTPCclusters(0x0);
	175
	176	Float_t chi2PerClusterITS = -1;
	177	if (nClustersITS!=0)
	178	chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
	179	Float_t chi2PerClusterTPC = -1;
	180	if (nClustersTPC!=0)
	181	chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
	182
	183	//
	184
	185	if(nClustersTPC < fMinNumberOfTPCClusters) continue;
	186	if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;
	187	if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;
	188	if(TMath::Abs(dcaZ) > fMaxDCAz) continue;
	189
	190	Short_t gCharge = tpcOnlyTrack->Charge();
	191	if(gCharge > 0) nPlus += 1;
	192	if(gCharge < 0) nMinus += 1;
	193	delete tpcOnlyTrack;
	194	} // track loop
	195
	196	if( nPlus == 0 && nMinus == 0) return;
	197
	198	// cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;
	199	fHistCentrality->Fill(nCentrality);
	200	fHistNPlusNMinus->Fill(nPlus,nMinus);
	201	fHistNMult->Fill(nPlus+nMinus);
	202
	203
	204	}//Vz cut
	205	}//Vy cut
206	}//Vx cut
207	}//Vz resolution
208	}//number of contributors
209	}//valid vertex
210	}//TPC analysis mode
211	}//centrality
212	//}//physics selection
213	}//ESD analysis level
214
215	//MC analysis
216	if(fAnalysisType.CompareTo("MC") == 0 ) {
217	AliMCEvent* mcEvent = MCEvent();
218	if (!mcEvent) {
219	Printf("ERROR: Could not retrieve MC event");
220	return;
221	}
222	AliStack* stack = mcEvent->Stack();
223	if (!stack) {
224	Printf("ERROR: Could not retrieve MC stack");
225	return;
226	}
227
228	fHistEventStats->Fill(1);
229	fHistEventStats->Fill(2);
230	fHistEventStats->Fill(3);
231	fHistEventStats->Fill(4); //analyzed events
232	for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
233	AliVParticle* track = mcEvent->GetTrack(iTracks);
234	if (!track) {
235	Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
236	continue;
237	}
238
239	if(!stack->IsPhysicalPrimary(iTracks)) continue;
240	if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
241	if(TMath::Abs(track->Eta()) > 0.8) continue;
242
243	Short_t gCharge = track->Charge();
244	if(gCharge > 0) nPlus += 1;
245	if(gCharge < 0) nMinus += 1;
246	}//particle loop
247	fHistNPlusNMinusMC->Fill(nPlus,nMinus);
248	fHistNMultMC->Fill(nPlus+nMinus);
249
250	}//MC analysis level
251
252
253	}
254
255	//________________________________________________________________________
256	void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {
257	// Draw result to the screen
258	// Called once at the end of the query
259
260	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
261	if (!fOutputList) {
262	printf("ERROR: Output list not available\n");
263	return;
264	}
265	}