[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEFluctuationAnalysisTask.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 "AliEbyEFluctuationAnalysisTask.h"

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

ClassImp(AliEbyEFluctuationAnalysisTask)

//________________________________________________________________________
AliEbyEFluctuationAnalysisTask::AliEbyEFluctuationAnalysisTask(const char *name) 
  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), 
    fHistEventStats(0), fHistCentrality(0),
    fHistNMultMC(0), fHistNPlusNMinusMC(0), 
    fESDtrackCuts(0),
    fAnalysisType("ESD"), fAnalysisMode("TPC"),
    fCentralityEstimator("V0M"), fCentralityBins20(kFALSE),
    fVxMax(3.), fVyMax(3.), fVzMax(10.) {
  // Constructor
  
  for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
    fHistNMult[iBin-1] = NULL;
    fHistNPlusNMinus[iBin-1] = NULL;
  }

  // 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 AliEbyEFluctuationAnalysisTask::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 == "ESD") {
    fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
			       nCentralityBins,0.5,nCentralityBins+0.5);
    fOutputList->Add(fHistCentrality);
    
    TString histName;
    for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
      histName = "fHistNMult"; histName += "Centrality"; histName += iBin; 
      fHistNMult[iBin-1] = new TH1F(histName.Data(), 
				    ";N_{mult.}",
				    500,0,3000);
      fOutputList->Add(fHistNMult[iBin-1]);
    }
    for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
      histName = "fHistNPlusNMinus"; histName += "Centrality"; histName += iBin; 
      fHistNPlusNMinus[iBin-1] = new TH2F(histName.Data(), 
					  ";N_{+};N_{-}",
					  2000,0.5,2000.5,2000,0.5,2000.5);  
      fOutputList->Add(fHistNPlusNMinus[iBin-1]);
    }
  }//ESD analysis
  else if(fAnalysisType == "MC") {
    TString histName = "fHistNMultMC";
    fHistNMultMC = new TH1F(histName.Data(), 
			    ";N_{mult.}",
			    600,0,6000);
    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
}

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

  // Post output data.
  //ESD analysis
  if(fAnalysisType == "ESD") {
    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
      Printf("Event accepted");

      //Centrality stuff
      AliCentrality *centrality = fESD->GetCentrality();
      Int_t nCentrality = 0;
      if(fCentralityBins20) 
	nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
      else
	nCentrality = centrality->GetCentralityClass10(fCentralityEstimator.Data());

      if((nCentrality < 0)||(nCentrality > 19)) return;
      
      if(fAnalysisMode == "TPC") {
	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->GetXv()) < fVxMax) {
		if(TMath::Abs(vertex->GetYv()) < fVyMax) {
		  if(TMath::Abs(vertex->GetZv()) < 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;
		    TObjArray *gTrackArray = 0;
		    if(fESDtrackCuts)
                      gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
                    if(gTrackArray) {
                      nAcceptedTracks = gTrackArray->GetEntries();

		      AliESDtrack* track = 0;
		      // Track loop to fill a pT spectrum
		      for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) {
			track = dynamic_cast<AliESDtrack *>(gTrackArray->At(iTracks));
			if (!track) {
			  printf("ERROR: Could not receive track %d\n", iTracks);
			  continue;
			}
			
			Short_t gCharge = track->Charge();
			if(gCharge > 0) nPlus += 1;
			if(gCharge < 0) nMinus += 1;
		      }//track loop
		    }//TObjArray valid object
		    //if((nCentrality >= 1)&&(nCentrality <= 20)) {
		    
		    fHistCentrality->Fill(nCentrality);
		    fHistNPlusNMinus[nCentrality-1]->Fill(nPlus,nMinus);
		    fHistNMult[nCentrality-1]->Fill(nPlus+nMinus);
		    //}
		  }//Vz cut
		}//Vy cut
	      }//Vx cut
	    }//Vz resolution
	  }//number of contributors
	}//valid vertex
      }//TPC analysis mode
      //}//physics selection
  }//ESD analysis level
  
  //MC analysis
  if(fAnalysisType == "MC") {
    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


  PostData(1, fOutputList);
}      

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