--- /dev/null
+/*************************************************************************
+* Copyright(c) 1998-2008, 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. *
+**************************************************************************/
+
+#include "Riostream.h" //needed as include
+#include "TChain.h"
+#include "TTree.h"
+#include "TFile.h" //needed as include
+#include "TList.h"
+
+
+class AliAnalysisTask;
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+
+#include "AliAnalysisTaskScalarProduct.h"
+#include "AliFlowEventSimpleMaker.h"
+#include "AliFlowAnalysisWithScalarProduct.h"
+
+// AliAnalysisTaskScalarProduct:
+//
+// analysis task for Scalar Product Method
+//
+// Author: Naomi van der Kolk (kolk@nikhef.nl)
+
+ClassImp(AliAnalysisTaskScalarProduct)
+
+//________________________________________________________________________
+AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct(const char *name) :
+ AliAnalysisTask(name, ""),
+ fESD(0),
+ fAOD(0),
+ fSP(0),
+ fEventMaker(0),
+ fAnalysisType("ESD")
+{
+ // Constructor
+ cout<<"AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct(const char *name)"<<endl;
+
+ // Define input and output slots here
+ // Input slot #0 works with a TChain
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TList container
+ DefineOutput(0, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScalarProduct::ConnectInputData(Option_t *)
+{
+ // Connect ESD or AOD here
+ // Called once
+ cout<<"AliAnalysisTaskScalarProduct::ConnectInputData(Option_t *)"<<endl;
+
+ TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+ if (!tree) {
+ Printf("ERROR: Could not read chain from input slot 0");
+ } else {
+ // Disable all branches and enable only the needed ones
+ if (fAnalysisType == "MC") {
+ cout<<"!!!!!reading MC kinematics only"<<endl;
+ // we want to process only MC
+ tree->SetBranchStatus("*", kFALSE);
+
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+ if (!esdH) {
+ Printf("ERROR: Could not get ESDInputHandler");
+ } else {
+ fESD = esdH->GetEvent();
+ }
+ }
+ else if (fAnalysisType == "ESD") {
+ cout<<"!!!!!reading the ESD only"<<endl;
+ tree->SetBranchStatus("*", kFALSE);
+ tree->SetBranchStatus("Tracks.*", kTRUE);
+
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+ if (!esdH) {
+ Printf("ERROR: Could not get ESDInputHandler");
+ } else
+ fESD = esdH->GetEvent();
+ }
+ else if (fAnalysisType == "AOD") {
+ cout<<"!!!!!reading the AOD only"<<endl;
+ AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+ if (!aodH) {
+ Printf("ERROR: Could not get AODInputHandler");
+ }
+ else {
+ fAOD = aodH->GetEvent();
+ }
+ }
+ else {
+ Printf("!!!!!Wrong analysis type: Only ESD, AOD and MC types are allowed!");
+ exit(1);
+
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScalarProduct::CreateOutputObjects()
+{
+ // Called once
+ cout<<"AliAnalysisTaskScalarProduct::CreateOutputObjects()"<<endl;
+
+ if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "MC")) {
+ cout<<"WRONG ANALYSIS TYPE! only ESD, AOD and MC are allowed."<<endl;
+ exit(1);
+ }
+
+ //event maker
+ fEventMaker = new AliFlowEventSimpleMaker();
+ //Analyser
+ fSP = new AliFlowAnalysisWithScalarProduct() ;
+
+ //output file
+ TString fName = "outputFromScalarProductAnalysis" ;
+ fName += fAnalysisType.Data() ;
+ fName += ".root" ;
+ fSP->SetHistFileName( fName.Data() );
+
+ fSP-> Init();
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScalarProduct::Exec(Option_t *)
+{
+ // Main loop
+ // Called for each event
+
+
+ // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
+ // This handler can return the current MC event
+
+ AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+ if (!eventHandler) {
+ Printf("ERROR: Could not retrieve MC event handler");
+ return;
+ }
+
+ AliMCEvent* mcEvent = eventHandler->MCEvent();
+ if (!mcEvent) {
+ Printf("ERROR: Could not retrieve MC event");
+ return;
+ }
+
+ Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
+
+ if (fAnalysisType == "MC") {
+ // analysis
+ AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent);
+ fSP->Make(fEvent);
+
+ delete fEvent;
+ }
+ else if (fAnalysisType == "ESD") {
+ if (!fESD) {
+ Printf("ERROR: fESD not available");
+ return;
+ }
+ Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
+
+ // analysis
+ AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
+ fSP->Make(fEvent);
+ delete fEvent;
+ }
+ else if (fAnalysisType == "AOD") {
+ if (!fAOD) {
+ Printf("ERROR: fAOD not available");
+ return;
+ }
+ Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
+
+ // analysis
+ AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
+ fSP->Make(fEvent);
+ delete fEvent;
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScalarProduct::Terminate(Option_t *)
+{
+ // Called once at the end of the query
+ fSP->Finish();
+ PostData(0,fSP->GetHistFile());
+
+ delete fSP;
+ delete fEventMaker;
+}
\ No newline at end of file
--- /dev/null
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id$ */
+
+#ifndef AliAnalysisTaskScalarProduct_H
+#define AliAnalysisTaskScalarProduct_H
+
+// AliAnalysisTaskScalarProduct:
+// analysis task for Scalar Product method
+// Author: Naomi van der Kolk (kolk@nikhef.nl)
+
+class AliESDEvent;
+class AliAODEvent;
+class AliFlowAnalysisWithScalarProduct;
+class AliFlowEventSimpleMaker;
+
+#include "TString.h"
+#include "AliAnalysisTask.h"
+
+class AliAnalysisTaskScalarProduct : public AliAnalysisTask {
+ public:
+ AliAnalysisTaskScalarProduct(const char *name = "AliAnalysisTaskScalarProduct");
+ virtual ~AliAnalysisTaskScalarProduct() {}
+
+ virtual void ConnectInputData(Option_t *);
+ virtual void CreateOutputObjects();
+ virtual void Exec(Option_t *option);
+ virtual void Terminate(Option_t *);
+
+ void SetAnalysisType(TString type) { this->fAnalysisType = type; }
+ TString GetAnalysisType() const { return this->fAnalysisType; }
+
+ private:
+ AliESDEvent *fESD; // ESD object
+ AliAODEvent *fAOD; // AOD object
+ AliFlowAnalysisWithScalarProduct* fSP; // analysis object
+ AliFlowEventSimpleMaker* fEventMaker; // FlowEventSimple maker object
+ TString fAnalysisType; // can be MC, ESD or AOD
+
+ ClassDef(AliAnalysisTaskScalarProduct, 1); // example of analysis
+};
+
+#endif
\ No newline at end of file
--- /dev/null
+/*************************************************************************
+* Copyright(c) 1998-2008, 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. *
+**************************************************************************/
+
+#define AliFlowAnalysisWithScalarProduct_cxx
+
+#include "Riostream.h" //needed as include
+#include "TFile.h" //needed as include
+#include "TMath.h"
+#include "TProfile.h"
+#include "TVector2.h"
+
+class TH1F;
+
+#include "AliFlowCommonConstants.h" //needed as include
+#include "AliFlowEventSimple.h"
+#include "AliFlowTrackSimple.h"
+#include "AliFlowCommonHist.h"
+//#include "AliFlowCommonHistResults.h"
+#include "AliFlowAnalysisWithScalarProduct.h"
+
+class AliFlowVector;
+
+// AliFlowAnalysisWithScalarProduct:
+// Description:
+// Maker to analyze Flow with the Scalar product method.
+//
+// author: N. van der Kolk (kolk@nikhef.nl)
+
+ClassImp(AliFlowAnalysisWithScalarProduct)
+
+ //-----------------------------------------------------------------------
+
+ AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
+ fEventNumber(0),
+ fEvent(0x0),
+ fTrack(0x0),
+ fDebug(kFALSE),
+ fHistFileName(0),
+ fHistFile(0),
+ fHistProUQ(0),
+ fCommonHists(0)
+{
+
+ // Constructor.
+ fQ.Set(0.,0.); // flow vector
+ fU.Set(0.,0.); // particle unit vector
+
+}
+ //-----------------------------------------------------------------------
+
+
+ AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
+ {
+ //destructor
+
+ }
+
+
+//-----------------------------------------------------------------------
+void AliFlowAnalysisWithScalarProduct::Init() {
+
+ //Define all histograms
+ cout<<"---Analysis with the Scalar Product Method---"<<endl;
+
+ Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
+ Double_t fPtMin = AliFlowCommonConstants::GetPtMin();
+ Double_t fPtMax = AliFlowCommonConstants::GetPtMax();
+
+ // analysis file (output)
+ fHistFile = new TFile(fHistFileName.Data(),"RECREATE") ;
+
+ fHistProUQ = new TProfile("Flow_UQ_SP","Flow_UQ_SP",fNbinsPt,fPtMin,fPtMax);
+ fHistProUQ->SetXTitle("p_t (GeV)");
+ fHistProUQ->SetYTitle("<uQ>");
+
+ fCommonHists = new AliFlowCommonHist("SP");
+ //fCommonHistsRes = new AliFlowCommonHistResults("SP");
+
+ fEventNumber = 0; //set number of events to zero
+
+}
+
+//-----------------------------------------------------------------------
+
+void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* fEvent) {
+
+ //Fill histogram
+ if (fEvent) {
+
+ //fill control histograms
+ fCommonHists->FillControlHistograms(fEvent);
+
+ //get the Q vector from the FlowEvent
+ fQ = fEvent->GetQ();
+ //Double_t fMult = fQ.GetMult();
+
+ //loop over the tracks of the event
+ Int_t fNumberOfTracks = fEvent->NumberOfTracks();
+ for (Int_t i=0;i<fNumberOfTracks;i++)
+ {
+
+ fTrack = fEvent->GetTrack(i) ;
+ if (fTrack){
+ if (fTrack->UseForDifferentialFlow()) {
+ Double_t fPhi = fTrack->Phi();
+
+ //calculate fU
+ Double_t fUX = TMath::Cos(2*fPhi);
+ Double_t fUY = TMath::Sin(2*fPhi);
+ fU.Set(fUX,fUY);
+ Double_t fModulus = fU.Mod();
+ if (fModulus!=0.) fU.Set(fUX/fModulus,fUY/fModulus); // make length 1
+ else cerr<<"fModulus is zero!"<<endl;
+
+ TVector2 fQm = fQ;
+ //subtrackt particle from the flowvector if used to define it
+ if (fTrack->UseForIntegratedFlow()) {
+ Double_t fQmX = fQm.X() - fUX;
+ Double_t fQmY = fQm.Y() - fUY;
+ fQm.Set(fQmX,fQmY);
+ }
+
+ //Double_t fUQ = scalar product of fU and fQm
+ Double_t fUQ = fU*fQm;
+ Double_t fPt = fTrack->Pt();
+ //fill the profile histogram
+ fHistProUQ->Fill(fPt,fUQ);
+ }
+ }//track selected
+ }//loop over tracks
+
+ fEventNumber++;
+ cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
+ }
+}
+
+ //--------------------------------------------------------------------
+void AliFlowAnalysisWithScalarProduct::Finish() {
+
+ //*************make histograms etc.
+ if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Terminate()"<<endl;
+
+ fHistProUQ->Draw();
+
+ // write to file
+ fHistFile->Write();
+
+ cout<<".....finished"<<endl;
+ }
+
+
\ No newline at end of file
--- /dev/null
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+
+#ifndef AliFlowAnalysisWithScalarProduct_H
+#define AliFlowAnalysisWithScalarProduct_H
+
+#include "TVector2.h" //called explicitly
+#include "AliFlowVector.h"
+#include "TString.h"
+
+class AliFlowTrackSimple;
+class AliFlowEventSimple;
+class AliFlowCommonHist;
+//class AliFlowCommonHistResults;
+
+class TProfile;
+class TFile;
+class Riostream;
+
+
+// Description: Maker to analyze Flow from the Scalar Product method.
+//
+// author: N. van der Kolk (kolk@nikhef.nl)
+
+
+class AliFlowAnalysisWithScalarProduct {
+
+ public:
+
+ AliFlowAnalysisWithScalarProduct(); //default constructor
+
+ virtual ~AliFlowAnalysisWithScalarProduct(); //destructor
+
+ void Init(); //defines variables and histograms
+ void Make(AliFlowEventSimple* fEvent); //calculates variables and fills histograms
+ void Finish(); //saves histograms
+
+ void SetDebug(Bool_t kt) { this->fDebug = kt ; }
+ Bool_t GetDebug() const { return this->fDebug ; }
+
+
+ // Output
+ void SetHistFileName(TString name) { this->fHistFileName = name ; } // Sets output file name
+ TString GetHistFileName() const { return this->fHistFileName ; } // Gets output file name
+ TFile* GetHistFile() const { return this->fHistFile ; } // Gets output file
+
+
+ private:
+ //cp const
+ //ass op
+
+ AliFlowVector fQ; // flow vector
+ TVector2 fU; // particle unit vector
+
+ Int_t fEventNumber; // event counter
+
+ AliFlowEventSimple* fEvent ; //!
+ AliFlowTrackSimple* fTrack ; //!
+
+ Bool_t fDebug ; //! flag for lyz analysis: more print statements
+
+ TString fHistFileName; //!
+ TFile* fHistFile; //!
+
+ TProfile* fHistProUQ; //!
+
+ AliFlowCommonHist* fCommonHists; //!
+ //AliFlowCommonHistResults* fCommonHistsRes; //!
+
+ ClassDef(AliFlowAnalysisWithScalarProduct,0) // macro for rootcint
+ };
+
+
+#endif
\ No newline at end of file
--- /dev/null
+// from CreateESDChain.C - instead of #include "CreateESDChain.C"
+TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 20, Int_t offset = 0) ;
+void LookupWrite(TChain* chain, const char* target) ;
+
+
+
+void runAliAnalysisTaskScalarProduct(Int_t nRuns = 2, TString type = "ESD", const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX/", Int_t offset = 0)
+
+{
+ TStopwatch timer;
+ timer.Start();
+
+ // include path (to find the .h files when compiling)
+ gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
+ gSystem->AddIncludePath("-I$ROOTSYS/include") ;
+
+ // load needed libraries
+ gSystem->Load("libTree.so");
+ gSystem->Load("libESD.so");
+ cerr<<"libESD loaded..."<<endl;
+ gSystem->Load("libANALYSIS.so");
+ cerr<<"libANALYSIS.so loaded..."<<endl;
+
+ gROOT->LoadMacro("AliFlowCommonConstants.cxx+");
+ gROOT->LoadMacro("AliFlowVector.cxx+");
+ gROOT->LoadMacro("AliFlowTrackSimple.cxx+");
+ gROOT->LoadMacro("AliFlowEventSimple.cxx+");
+ gROOT->LoadMacro("AliFlowEventSimpleMaker.cxx+");
+ gROOT->LoadMacro("AliFlowCommonHist.cxx+");
+ //gROOT->LoadMacro("AliFlowCommonHistResults.cxx+");
+ gROOT->LoadMacro("AliFlowAnalysisWithScalarProduct.cxx+");
+ gROOT->LoadMacro("AliAnalysisTaskScalarProduct.cxx+");
+
+ // create the TChain. CreateESDChain() is defined in CreateESDChain.C
+ TChain* chain = CreateESDChain(dataDir, nRuns, offset);
+ cout<<"chain ("<<chain<<")"<<endl;
+
+ //____________________________________________//
+ // Make the analysis manager
+ AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+
+ if (type == "ESD"){
+ AliVEventHandler* esdH = new AliESDInputHandler;
+ mgr->SetInputEventHandler(esdH); }
+
+ if (type == "AOD"){
+ AliVEventHandler* aodH = new AliAODInputHandler;
+ mgr->SetInputEventHandler(aodH); }
+
+ AliMCEventHandler *mc = new AliMCEventHandler();
+ mgr->SetMCtruthEventHandler(mc);
+ //____________________________________________//
+ // 1st MC EP task
+ AliAnalysisTaskScalarProduct *task1 = new AliAnalysisTaskScalarProduct("TaskScalarProduct");
+ task1->SetAnalysisType(type);
+ mgr->AddTask(task1);
+
+ // Create containers for input/output
+ AliAnalysisDataContainer *cinput1 =
+ mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
+ AliAnalysisDataContainer *coutput1 =
+ mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer);
+
+ //____________________________________________//
+ mgr->ConnectInput(task1,0,cinput1);
+ mgr->ConnectOutput(task1,0,coutput1);
+
+ if (!mgr->InitAnalysis()) return;
+ mgr->PrintStatus();
+ mgr->StartAnalysis("local",chain);
+
+ timer.Stop();
+ timer.Print();
+}
+
+
+// Helper macros for creating chains
+// from: CreateESDChain.C,v 1.10 jgrosseo Exp
+
+TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+ // creates chain of files in a given directory or file containing a list.
+ // In case of directory the structure is expected as:
+ // <aDataDir>/<dir0>/AliESDs.root
+ // <aDataDir>/<dir1>/AliESDs.root
+ // ...
+
+ if (!aDataDir)
+ return 0;
+
+ Long_t id, size, flags, modtime;
+ if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+ {
+ printf("%s not found.\n", aDataDir);
+ return 0;
+ }
+
+ TChain* chain = new TChain("esdTree");
+ TChain* chaingAlice = 0;
+
+ if (flags & 2)
+ {
+ TString execDir(gSystem->pwd());
+ TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+ TList* dirList = baseDir->GetListOfFiles();
+ Int_t nDirs = dirList->GetEntries();
+ gSystem->cd(execDir);
+
+ Int_t count = 0;
+
+ for (Int_t iDir=0; iDir<nDirs; ++iDir)
+ {
+ TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+ if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+ continue;
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ TString presentDirName(aDataDir);
+ presentDirName += "/";
+ presentDirName += presentDir->GetName();
+ chain->Add(presentDirName + "/AliESDs.root/esdTree");
+ cerr<<presentDirName<<endl;
+ }
+
+ }
+ else
+ {
+ // Open the input stream
+ ifstream in;
+ in.open(aDataDir);
+
+ Int_t count = 0;
+
+ // Read the input list of files and add them to the chain
+ TString esdfile;
+ while(in.good()) {
+ in >> esdfile;
+ if (!esdfile.Contains("root")) continue; // protection
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ // add esd file
+ chain->Add(esdfile);
+ }
+
+ in.close();
+ }
+
+ return chain;
+}
+
+void LookupWrite(TChain* chain, const char* target)
+{
+ // looks up the chain and writes the remaining files to the text file target
+
+ chain->Lookup();
+
+ TObjArray* list = chain->GetListOfFiles();
+ TIterator* iter = list->MakeIterator();
+ TObject* obj = 0;
+
+ ofstream outfile;
+ outfile.open(target);
+
+ while ((obj = iter->Next()))
+ outfile << obj->GetTitle() << "#AliESDs.root" << endl;
+
+ outfile.close();
+
+ delete iter;
+}