From: mfloris Date: Fri, 5 Nov 2010 09:07:04 +0000 (+0000) Subject: Main Updates: X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=e0376287fa7c2c6c4eaf2ea72b205850f7c8b7cf Main Updates: - AliAnalysisMultPbCentralitySelector: when running on data, uses AliESDcentrality by default, when running on MC a cut on the selected tracks - AliAnalysisMultPbTrackHistoManager&AliAnalysisTaskMultPbTracks.cxx: added histograms for particle species, multiplicity histogram, changed pt binning (same as 900 GeV pp now) - AliAnalysisTaskTriggerStudy: added venn-like histo - AddTaskMultPbPbTracks.C: forgot in the initial commit - run.C&run.sh: changes to run on grid (not yet working), loading extra classes through a TList, possibility to give a custom prefix to the output folder. - correct.C: minor chages for specific checks --- diff --git a/PWG0/multPbPb/AddTaskMultPbPbTracks.C b/PWG0/multPbPb/AddTaskMultPbPbTracks.C new file mode 100644 index 00000000000..a99b2b5ce65 --- /dev/null +++ b/PWG0/multPbPb/AddTaskMultPbPbTracks.C @@ -0,0 +1,66 @@ +AliAnalysisTaskMultPbTracks * AddTaskMultPbPbTracks(const char * outfilename, AliESDtrackCuts * esdTrackCuts = 0) +{ + // TODO: add some parameters to set the centrality for this task, and maybe the name of the task + // TODO: shall I use the same file and different dirs for the different centralities? + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskPhysicsSelection", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskPhysicsSelection", "This task requires an input event handler"); + return NULL; + } + TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + + if (inputDataType != "ESD") { + Printf("ERROR! This task can only run on ESDs!"); + } + + // Configure analysis + //=========================================================================== + + + + AliAnalysisTaskMultPbTracks *task = new AliAnalysisTaskMultPbTracks("TaskMultPbTracks"); + mgr->AddTask(task); + + // Set Cuts + if (!esdTrackCuts) + { + printf("ERROR: esdTrackCuts could not be created\n"); + return; + } + task->SetTrackCuts(esdTrackCuts); + + // TODO: + // IO into folders in a file? + // FIXME: + // Add centrality selection configuration (included arguments) + + // Set I/O + AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("cmultPbTracksOutHM", + AliAnalysisMultPbTrackHistoManager::Class(), + AliAnalysisManager::kOutputContainer, + outfilename); + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("cmultPbTracksOutCT", + AliESDtrackCuts::Class(), + AliAnalysisManager::kOutputContainer, + outfilename); + // AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cmultPbTracksOutCM", + // AliAnalysisMultPbCentralitySelector::Class(), + // AliAnalysisManager::kOutputContainer, + // outfilename); + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task,1,coutput1); + mgr->ConnectOutput(task,2,coutput2); + // mgr->ConnectOutput(task,3,coutput3); + + return task; +} diff --git a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx index 9a3bd65b1a4..6029b87a23b 100644 --- a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx +++ b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx @@ -1,4 +1,45 @@ #include "AliAnalysisMultPbCentralitySelector.h" +#include "AliESDtrackCuts.h" +#include "AliESDCentrality.h" +#include "AliESDEvent.h" +#include "AliLog.h" ClassImp(AliAnalysisMultPbCentralitySelector) +Bool_t AliAnalysisMultPbCentralitySelector::IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts) { + + // Centrality selection + // On MC cuts on the number of good tracks, + // On data cuts using AliESDCentrality and the cut requested in ntracks + if(fIsMC) { + if(!trackCuts){ + AliFatal("Track cuts object is invalid"); + } + if (trackCuts->CountAcceptedTracks(aEsd) < fMultMin) return kFALSE; + if (trackCuts->CountAcceptedTracks(aEsd) > fMultMax) return kFALSE; + } + + AliESDCentrality *centrality = aEsd->GetCentrality(); + if(!centrality) { + AliFatal("Centrality object not available"); + } + else { + Int_t centrBin = centrality->GetCentralityClass5(fCentrEstimator.Data()) ; + if (centrBin != fCentrBin && fCentrBin != -1) return kFALSE; + } + + return kTRUE; + +} + +void AliAnalysisMultPbCentralitySelector::Print(Option_t* option ) const { + // Print some information + + Printf("AliAnalysisMultPbCentralitySelector"); + Printf(" - Centrality estimator [%s]",fCentrEstimator.Data()); + + if ( fIsMC ) { + Printf("Running on Monte Carlo, actual cut was on tracks multiplicity [%d - %d]",fMultMin,fMultMax); + } + +} diff --git a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h index eb07d7cc37f..0053722fd7e 100644 --- a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h +++ b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h @@ -24,30 +24,34 @@ class TH1F; class TCollection; class AliTriggerAnalysis; class AliAnalysisTaskSE; - +class AliESDtrackCuts; class AliAnalysisMultPbCentralitySelector : public AliAnalysisCuts { public: - AliAnalysisMultPbCentralitySelector() : fIsMC (0) {;} + AliAnalysisMultPbCentralitySelector() : fIsMC (0), fCentrEstimator(""), fCentrBin(-1), fMultMin(0), fMultMax(1000000) {;} virtual ~AliAnalysisMultPbCentralitySelector(){} // AliAnalysisCuts interface - virtual UInt_t GetSelectionMask(const TObject* obj) { return IsCentralityBinSelected((const AliESDEvent*) obj); } + virtual UInt_t GetSelectionMask(const TObject* obj) { return (UInt_t) IsCentralityBinSelected((AliESDEvent*) obj, NULL); } virtual Bool_t IsSelected(TList*) { AliFatal("Not implemented"); return kFALSE; } - virtual Bool_t IsSelected(TObject* obj) {return IsCentralityBinSelected ( (AliESDEvent*) obj);} - - UInt_t IsCentralityBinSelected(const AliESDEvent* aEsd){ return kTRUE;} + virtual Bool_t IsSelected(TObject* obj) {return (UInt_t) IsCentralityBinSelected ( (AliESDEvent*) obj, NULL);} - void SetAnalyzeMC(Bool_t flag = kTRUE) { fIsMC = flag; } + Bool_t IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts); - virtual void Print(Option_t* option = "") const { Printf ("Multiplitity Selector [AliAnalysisMultPbCentralitySelector] [%s]", option);} + void SetAnalyzeMC(Bool_t flag = kTRUE, Double_t multMin = 0, Double_t multMax=10000) { fIsMC = flag; fMultMin = multMin; fMultMax = multMax; } + void SetCentralityEstimator(const char * estimator) { fCentrEstimator = estimator; } + + virtual void Print(Option_t* option = "") const ; virtual Long64_t Merge(TCollection* list){list->GetEntries();return 0;} protected: Bool_t fIsMC; // flag if MC is analyzed - + TString fCentrEstimator; // Centrality estimator for AliESDCentrality + TString fCentrBin; // centrality bin to be selected + Int_t fMultMin ; // Minimum multiplicity, because on MC we cut on tracks rather than on the estimator + Int_t fMultMax ; // Maximum multiplicity, because on MC we cut on tracks rather than on the estimator ClassDef(AliAnalysisMultPbCentralitySelector, 1) private: diff --git a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx index dda337c4572..ca2aa7a3388 100644 --- a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx +++ b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx @@ -4,6 +4,8 @@ #include "TH3D.h" #include "TH1I.h" #include "TROOT.h" +#include "TMCProcess.h" + #include using namespace std; ClassImp(AliAnalysisMultPbTrackHistoManager) @@ -12,6 +14,7 @@ const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[] = { "All E const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim", "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"}; const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[] = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"}; +const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake"}; @@ -58,6 +61,20 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) { } +TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) { + // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it. + + TString name = kHistoPrefix[id]; + name += "Mult"; + TH1D * h = (TH1D*) GetHisto(name.Data()); + if (!h) { + h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id])); + } + + return h; + +} + TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, Float_t minEta, Float_t maxEta, Float_t minVz, Float_t maxVz, @@ -84,7 +101,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, if (gROOT->FindObjectAny(hname.Data())){ - AliError(Form("An object called %s already exists",hname.Data())); + AliError(Form("An object called %s already exists,adding suffix",hname.Data())); + hname += "_2"; } TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E"); @@ -119,7 +137,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id, hname = hname + "_Vz_" + long (min1) + "_" + long(max1) + "_" + long (min2) + "_" + long(max2); if (gROOT->FindObjectAny(hname.Data())){ - AliError(Form("An object called %s already exists",hname.Data())); + AliError(Form("An object called %s already exists, adding suffix",hname.Data())); + hname+="_2"; } TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E"); @@ -153,7 +172,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id, hname = hname + "_Eta_" + long (min1) + "_" + long(max1) + "_" + long (min2) + "_" + long(max2); if (gROOT->FindObjectAny(hname.Data())){ - AliError(Form("An object called %s already exists",hname.Data())); + AliError(Form("An object called %s already exists, adding suffix",hname.Data())); + hname+="_2"; } TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E"); @@ -161,6 +181,36 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id, return h; } +TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) { + + // Returns histogram with particle specties + + TString name = TString(kHistoPrefix[id])+"_Species"; + + TH1D * h = (TH1D*) GetHisto(name); + if (!h) { + name+=fHNameSuffix; + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + AliInfo(Form("Booking histo %s",name.Data())); + + h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5); + Int_t nbin = kPNoProcess+1; + for(Int_t ibin = 0; ibin < nbin; ibin++){ + h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]); + } + TH1::AddDirectory(oldStatus); + fList->Add(h); + + + } + return h; + + +} + + TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() { // Returns histogram with event statistiscs (processed events at each step) @@ -189,7 +239,7 @@ void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * op if (!h->InheritsFrom("TH1")) { AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName())); } - cout << "Scaling " << h->GetName() << " " << nev << endl; + AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev)); h->Scale(1./nev,option); } @@ -208,15 +258,39 @@ TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, c AliInfo(Form("Booking %s",hname.Data())); + // Binning from Jacek task + const Int_t nptbins = 68; + const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0}; + + const Int_t netabins=20; + Double_t binsEta[netabins+1]; + Float_t minEta = -1; + Float_t maxEta = 1; + Float_t etaStep = (maxEta-minEta)/netabins; + for(Int_t ibin = 0; ibin < netabins; ibin++){ + binsEta[ibin] = minEta + ibin*etaStep; + binsEta[ibin+1] = minEta + ibin*etaStep + etaStep; + } + + const Int_t nvzbins=10; + Double_t binsVz[nvzbins+1]; + Float_t minVz = -10; + Float_t maxVz = 10; + Float_t vzStep = (maxVz-minVz)/nvzbins; + for(Int_t ibin = 0; ibin < nvzbins; ibin++){ + binsVz[ibin] = minVz + ibin*vzStep; + binsVz[ibin+1] = minVz + ibin*vzStep + vzStep; + } + TH3D * h = new TH3D (hname,title, - 50,0,10, // Pt - 20,-1, 1, // Eta - 10,-10,10 // Vz + nptbins, binsPt, + netabins, binsEta, + nvzbins, binsVz ); - h->SetYTitle("#eta"); h->SetXTitle("p_{T}"); + h->SetYTitle("#eta"); h->SetZTitle("V_{z} (cm)"); h->Sumw2(); @@ -248,6 +322,28 @@ TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const TH1::AddDirectory(oldStatus); return h; } +TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) { + // Books a multiplicity histo + + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + TString hname = name; + hname+=fHNameSuffix; + + AliInfo(Form("Booking %s",hname.Data())); + + + TH1D * h = new TH1D (hname,title, 600, 0,6000); + + h->SetXTitle("N tracks"); + h->Sumw2(); + + fList->Add(h); + + TH1::AddDirectory(oldStatus); + return h; +} TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() { // Books histogram with event statistiscs (processed events at each step) diff --git a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h index 4e74ce6ae5d..6f0591bcaa7 100644 --- a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h +++ b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h @@ -42,6 +42,9 @@ public: TH1I * GetHistoStats(); TH1D * GetHistoDCA(Histo_t id); + TH1D * GetHistoMult(Histo_t id); + + TH1D * GetHistoSpecies(Histo_t id); // Misch utils void ScaleHistos (Double_t nev, Option_t * option=""); @@ -52,6 +55,7 @@ public: TH3D * BookHistoPtEtaVz(const char * name, const char * title); TH1D * BookHistoDCA(const char * name, const char * title); TH1I * BookHistoStats(); + TH1D * BookHistoMult(const char * name, const char * title); // TH1 * GetHisto(const char * name); @@ -61,6 +65,7 @@ private: static const char * kStatStepNames[]; // names of the step hist static const char * kHistoPtEtaVzNames[]; // names of the 3D histograms pt/eta/vz static const char * kHistoDCANames[]; // names of the DCA histograms + static const char * kHistoPrefix[]; // prefix for histo names // FIXME: remove the others and keep only this TString fHNameSuffix; // Suffix added to all histo names. Useful if you have in the same session e.g. MC and data. AliAnalysisMultPbTrackHistoManager& operator=(const AliAnalysisMultPbTrackHistoManager& task); diff --git a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx index 345417bc7ca..11532c803c3 100644 --- a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx +++ b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx @@ -3,7 +3,8 @@ // Author: Michele Floris, CERN // TODO: // - Add chi2/cluster plot for primary, secondaries and fakes - +// FIXME: +// - re-implement centrality estimator to cut on tracks and multiplicity #include "AliAnalysisTaskMultPbTracks.h" #include "AliESDInputHandler.h" @@ -32,7 +33,7 @@ AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks() DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class()); DefineOutput(2, AliESDtrackCuts::Class()); - // DefineOutput(2, TH1I::Class()); + // DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class()); fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis"); if(fIsMC) fHistoManager->SetSuffix("MC"); @@ -48,6 +49,7 @@ AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name) DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class()); DefineOutput(2, AliESDtrackCuts::Class()); + // DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class()); fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis"); if(fIsMC) fHistoManager->SetSuffix("MC"); @@ -102,10 +104,12 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) /* PostData(0) is taken care of by AliAnalysisTaskSE */ PostData(1,fHistoManager); PostData(2,fTrackCuts); + // PostData(3,fCentralityEstimator); // Cache histogram pointers - static TH3D * hTracks[AliAnalysisMultPbTrackHistoManager::kNHistos]; - static TH1D * hDCA [AliAnalysisMultPbTrackHistoManager::kNHistos]; + static TH3D * hTracks [AliAnalysisMultPbTrackHistoManager::kNHistos]; + static TH1D * hDCA [AliAnalysisMultPbTrackHistoManager::kNHistos]; + static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos]; hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen ); hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec ); hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ); @@ -120,6 +124,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat ); hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak); + hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen ); fESD = dynamic_cast(fInputEvent); if (strcmp(fESD->ClassName(),"AliESDEvent")) { @@ -130,13 +135,15 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll); // Centrality selection + // Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD); + // if(!isCentralitySelected) return; + AliESDCentrality *centrality = fESD->GetCentrality(); if(!centrality) { - AliError("Centrality object not available"); // FIXME AliFatal here after debug + AliFatal("Centrality object not available"); } else { Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ; - cout << fCentralityEstimator.Data() << " BIN: " << centrBin << endl; if (centrBin != fCentrBin && fCentrBin != -1) return; } @@ -150,27 +157,25 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) AliError("No MC info found"); } else { - //loop on the MC event - // Int_t nMCTracks = fMCEvent->GetNumberOfTracks(); - Int_t offset = fMCEvent->GetPrimaryOffset(); - Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset; - for (Int_t ipart=offset; ipartGetNumberOfPrimaries(); + Int_t nPhysicalPrimaries = 0; + for (Int_t ipart=0; ipartGetTrack(ipart); // We don't care about neutrals and non-physical primaries if(mcPart->Charge() == 0) continue; - // FIXME: add kTransportBit (uncomment below) - if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue; - //check if current particle is a physical primary - // Bool_t physprim=fMCEvent->IsPhysicalPrimary(label); - // if (!physprim) continue; - // if (!track) return kFALSE; - // Bool_t transported = mcPart->Particle()->TestBit(kTransportBit); - // if(!transported) return kFALSE; + if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue; + nPhysicalPrimaries++; + // Fill species histo + fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID()); + + // Get MC vertex //FIXME: which vertex do I take for MC? TArrayF vertex; @@ -181,9 +186,12 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv); } + hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries); + } } + // FIXME: shall I take the primary vertex? const AliESDVertex* vtxESD = fESD->GetPrimaryVertex(); if(!vtxESD) return; @@ -242,7 +250,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) if (fIsMC) { // Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!! Int_t label = esdTrack->GetLabel(); // - AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(label); + AliMCParticle *mcPart = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label); if (!mcPart) { if(accepted) hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ()); @@ -250,8 +258,9 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA); } else { - if(fMCEvent->Stack()->IsPhysicalPrimary(label)) { - // FIXME add kTransportBit + if(IsPhysicalPrimaryAndTransportBit(label)) { + // Fill species histo + fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID()); if(accepted) hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ()); if(acceptedNoDCA) @@ -266,11 +275,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *) mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth)))); } if(mfl==3){ // strangeness + fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID()); if(accepted) hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ()); if(acceptedNoDCA) hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA); }else{ // material + fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID()); if(accepted) hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ()); if(acceptedNoDCA) @@ -295,3 +306,13 @@ void AliAnalysisTaskMultPbTracks::Terminate(Option_t *){ } +Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) { + + Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart); + if (!physprim) return kFALSE; + Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit); + if(!transported) return kFALSE; + + return kTRUE; + +} diff --git a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h index 8999008032e..3bd797949fe 100644 --- a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h +++ b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h @@ -29,12 +29,14 @@ public: AliAnalysisTaskMultPbTracks(const char * name); AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) ; ~AliAnalysisTaskMultPbTracks(); + //void SetCentralitySelector(AliAnalysisMultPbCentralitySelector * centr) { fCentrSelector=centr;} void SetTrackCuts(AliESDtrackCuts * cuts) { fTrackCuts = cuts;} void SetCentralityBin(Int_t bin = 0) { fCentrBin = bin; } void SetCentralityEstimator(const char * centr) { fCentralityEstimator = centr; } void SetIsMC(Bool_t flag=kTRUE) { fIsMC = flag;} AliAnalysisMultPbTrackHistoManager * GetHistoManager() { return fHistoManager;} + Bool_t IsPhysicalPrimaryAndTransportBit(Int_t ipart) ; virtual void UserCreateOutputObjects(); virtual void UserExec(Option_t *option); @@ -48,6 +50,7 @@ private: AliESDEvent * fESD; //! ESD object AliVEvent* fEvent; // TList * fListHisto; // list of output object AliAnalysisMultPbTrackHistoManager * fHistoManager; // wrapper for the list, takes care of merging + histo booking and getters + // AliAnalysisMultPbCentralitySelector * fCentrSelector; // centrality selector Int_t fCentrBin; // centrality bin selected (5% XS percentiles) TString fCentralityEstimator; // Name of the centrality estimator, for AliESDCentrality AliESDtrackCuts * fTrackCuts; // track cuts diff --git a/PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx b/PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx index 4b88d3a15cf..ad5f5fb8b16 100644 --- a/PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx +++ b/PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx @@ -94,10 +94,6 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *) AliFatal("Not processing ESDs"); } - // FIXME: two options here: either we add a loop setting the name of - // the histos with the trigger class in them (there may be a global - // SetHistoNamePrefix) or I put a cut to select only collision - // classes // get the multiplicity object const AliMultiplicity* mult = fESD->GetMultiplicity(); @@ -209,19 +205,9 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *) // // We don't care about neutrals and non-physical primaries // if(mcPart->Charge() == 0) continue; - // // FIXME: add kTransportBit (uncomment below) - // if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue; - - // //check if current particle is a physical primary - // // Bool_t physprim=fMCEvent->IsPhysicalPrimary(label); - // // if (!physprim) continue; - // // if (!track) return kFALSE; - // // Bool_t transported = mcPart->Particle()->TestBit(kTransportBit); - // // if(!transported) return kFALSE; - + // PHYSICAL PRIMARY // // Get MC vertex - // //FIXME: which vertex do I take for MC? - // TArrayF vertex; + // TArrayF vertex; // fMCEvent->GenEventHeader()->PrimaryVertex(vertex); // Float_t zv = vertex[2]; // // Float_t zv = vtxESD->GetZ(); @@ -313,19 +299,19 @@ void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const Bool_t fo2 = nFastOrOffline>=2; - if(fo2 && v0A && v0C && OM2) {cout << "Bin5: " << ibin << endl; h->Fill(1);} - if(fo2 && !v0A && v0C && OM2) {cout << "Bin6: " << ibin << endl; h->Fill(2);} - if(fo2 && v0A && !v0C && OM2) {cout << "Bin7: " << ibin << endl; h->Fill(3);} - if(fo2 && v0A && v0C && !OM2) {cout << "Bin8: " << ibin << endl; h->Fill(4);} - if(fo2 && v0A && !v0C && !OM2) {cout << "Bin9: " << ibin << endl; h->Fill(5);} - if(fo2 && !v0A && v0C && !OM2) {cout << "Bin10: " << ibin << endl; h->Fill(6);} - if(fo2 && !v0A && !v0C && OM2) {cout << "Bin11: " << ibin << endl; h->Fill(7);} - if(!fo2 && v0A && !v0C && OM2) {cout << "Bin12: " << ibin << endl; h->Fill(8);} - if(!fo2 && !v0A && v0C && OM2) {cout << "Bin13: " << ibin << endl; h->Fill(9);} - if(!fo2 && v0A && v0C && !OM2) {cout << "Bin14: " << ibin << endl; h->Fill(10);} - if(fo2 && !v0A && !v0C && !OM2) {cout << "Bin1: " << ibin << endl; h->Fill(11);} - if(!fo2 && v0A && !v0C && !OM2) {cout << "Bin2: " << ibin << endl; h->Fill(12);} - if(!fo2 && !v0A && v0C && !OM2) {cout << "Bin3: " << ibin << endl; h->Fill(13);} - if(!fo2 && !v0A && !v0C && OM2) {cout << "Bin4: " << ibin << endl; h->Fill(14);} + if(fo2 && v0A && v0C && OM2) { h->Fill(1);} + if(fo2 && !v0A && v0C && OM2) { h->Fill(2);} + if(fo2 && v0A && !v0C && OM2) { h->Fill(3);} + if(fo2 && v0A && v0C && !OM2) { h->Fill(4);} + if(fo2 && v0A && !v0C && !OM2) { h->Fill(5);} + if(fo2 && !v0A && v0C && !OM2) { h->Fill(6);} + if(fo2 && !v0A && !v0C && OM2) { h->Fill(7);} + if(!fo2 && v0A && !v0C && OM2) { h->Fill(8);} + if(!fo2 && !v0A && v0C && OM2) { h->Fill(9);} + if(!fo2 && v0A && v0C && !OM2) { h->Fill(10);} + if(fo2 && !v0A && !v0C && !OM2) { h->Fill(11);} + if(!fo2 && v0A && !v0C && !OM2) { h->Fill(12);} + if(!fo2 && !v0A && v0C && !OM2) { h->Fill(13);} + if(!fo2 && !v0A && !v0C && OM2) { h->Fill(14);} } diff --git a/PWG0/multPbPb/correct.C b/PWG0/multPbPb/correct.C index 532dda7863e..244dfab3bf7 100644 --- a/PWG0/multPbPb/correct.C +++ b/PWG0/multPbPb/correct.C @@ -24,50 +24,57 @@ TH1D * gHistoCompoments[kHistoFitCompoments]; void LoadLibs(); void LoadData(TString dataFolder, TString correctionFolder); void SetStyle(); -Double_t CheckSecondaries(); +void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial); +void CheckVz(); void ShowAcceptanceInVzSlices() ; TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ; - TH1D * GetCumulativeHisto (TH1 * h) ; +TH1D * GetCumulativeHisto (TH1 * h) ; static Double_t HistoSum(const double * x, const double* p); TF1 * GetFunctionHistoSum() ; TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ; TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ; TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ; +void PrintCanvas(TCanvas* c,const TString formats) ; +// global switches +Bool_t doPrint=kFALSE;// disable PrintCanvas - -void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correctionFolder = "output/LHC10a8_104867/") { +void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TString correctionFolder = "./output/LHC10g2a_130844_V0M_bin_10/") { // Load stuff and set some styles LoadLibs(); LoadData(dataFolder,correctionFolder); SetStyle(); - // ShowAcceptanceInVzSlices(); - // return; + ShowAcceptanceInVzSlices(); + return; // TODO add some cool printout for cuts and centrality selection - Double_t fractionWeak = CheckSecondaries(); - cout << "Rescaling weak correction: " << fractionWeak << endl; + CheckVz(); + + Double_t fractionWeak = 1, fractionMaterial=1; + // CheckSecondaries(fractionWeak, fractionMaterial); + cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10)->Clone("hDataPt"); - TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,-10,10); - TH1D * hMCPtRec = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10); - TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, -0.5,0.5,-10,10); - TH1D * hMCPtSeM = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, -0.5,0.5,-10,10); - TH1D * hMCPtSeW = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, -0.5,0.5,-10,10); - TH1D * hMCPtFak = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, -0.5,0.5,-10,10); + TH1D * hDataPt = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10)->Clone("hDataPt"); + TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,-10,10); + TH1D * hMCPtRec = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10); + TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, -0.5,0.5,-10,10); + TH1D * hMCPtSeM = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, -0.5,0.5,-10,10); + TH1D * hMCPtSeW = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, -0.5,0.5,-10,10); + TH1D * hMCPtFak = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, -0.5,0.5,-10,10); TCanvas * cdata = new TCanvas ("cData", "Data"); - + cdata->SetLogy(); hDataPt->Draw(); // hMCPtRec->Draw("same"); TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo"); + cMC->SetLogy(); cMC->cd(); hMCPtGen ->Draw(); hMCPtRec ->Draw("same"); @@ -90,6 +97,7 @@ void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correct TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hEffPt"); hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B"); + hCorSeM->Scale(fractionMaterial);// rescale material correction hCorSeM->Multiply(hDataPt); TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hEffPt"); @@ -125,7 +133,7 @@ void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correct hDataGen->Draw("same"); } -Double_t CheckSecondaries() { +void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) { // Returns the fraction you need to rescale the secondaries from weak decays for. // Some shorthands @@ -149,35 +157,54 @@ Double_t CheckSecondaries() { TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions"); c1->SetLogy(); // Draw all - hDataDCA->Draw(); - hMCDCARec ->Draw("same"); - hMCDCAPri ->Draw("same"); - hMCDCASW ->Draw("same"); - hMCDCASM ->Draw("same"); - hMCDCAFak ->Draw("same"); - return 1; + // hDataDCA->Draw(); + // hMCDCARec ->Draw("same"); + // hMCDCAPri ->Draw("same"); + // hMCDCASW ->Draw("same"); + // hMCDCASM ->Draw("same"); + // hMCDCAFak ->Draw("same"); + // return; + // Fit the DCA distribution. Uses a TF1 made by summing histograms TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak"); // hMCPrimSMFak->Add(hMCDCASM); hMCPrimSMFak->Add(hMCDCAFak); + // Set the components which are used in HistoSum, the static + // function for GetFunctionHistoSum gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone(); gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone(); gHistoCompoments[2] = (TH1D*) hMCDCASM->Clone(); TF1 * fHistos = GetFunctionHistoSum(); - fHistos->SetParameters(1,1); + fHistos->SetParameters(1,1,1); + fHistos->SetLineColor(kRed); + // Fit! hDataDCA->Fit(fHistos,"","",0,200); + // Rescale the components and draw to see how it looks like hMCPrimSMFak->Scale(fHistos->GetParameter(0)); hMCDCASW ->Scale(fHistos->GetParameter(1)); hMCDCASM ->Scale(fHistos->GetParameter(2)); hMCPrimSMFak->Draw("same"); hMCDCASW ->Draw("same"); hMCDCASM ->Draw("same"); - return fHistos->GetParameter(1)/fHistos->GetParameter(0); + // compute scaling factors + fracWeak = hMCDCASW->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral()); + fracMaterial = hMCDCASM->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral()); } +void CheckVz() { + // compares the Vz distribution in data and in MC + TCanvas * c1 = new TCanvas("cVz", "Vertex Z distribution"); + c1->cd(); + TH1D * hData = hManData->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec ); + TH1D * hCorr = hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec ); + hCorr->Draw(""); + hData->Draw("same"); + +} + void LoadLibs() { gSystem->Load("libVMC"); @@ -297,8 +324,10 @@ void LoadData(TString dataFolder, TString correctionFolder){ Int_t irowGoodTrigger = 1; if (hEvStatCorr && hEvStatData) { // hManData->ScaleHistos(75351.36/1.015);// Nint for run 104892 estimated correcting for the trigger efficiency, multiplied for the physics selection efficiency which I'm not correcting for the time being - hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger)); - hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger)); + // hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger)); + // hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger)); + hManData->ScaleHistos(hManData->GetHistoStats()->GetBinContent(2)); + hManCorr->ScaleHistos(hManCorr->GetHistoStats()->GetBinContent(2)); } else { cout << "WARNING!!! ARBITRARY SCALING" << endl; hManData->ScaleHistos(1000); @@ -384,22 +413,32 @@ TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) { void ShowAcceptanceInVzSlices() { TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz"); - for(Int_t ivz = -10; ivz < 10; ivz+=4){ - TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , -0.5,0.5,ivz,ivz+4); - TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,ivz,ivz+4); + for(Int_t ivz = -10; ivz < -6; ivz+=2){ + ivz=0;//FIXME + Bool_t first = kTRUE; + TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , -0.5,0.5,ivz,ivz+2); + TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,ivz,ivz+2); + TH1D * hMCPtPriD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , -0.5,0.5,ivz,ivz+2); + TH1D * hMCPtGenD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,ivz,ivz+2); // hEff= hMCPtGen; - TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+4)); + TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+2)); hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B"); + TH1D * hEffD = (TH1D*)hMCPtPriD->Clone(Form("hEffD_vz_%d_%d",ivz,ivz+2)); + hEffD->Divide(hMCPtPriD,hMCPtGenD,1,1,"B"); + hEffD->SetLineColor(kRed); cout << "ivz " << ivz << endl; - if(ivz < -9) { + if(first) { + first=kFALSE; cout << "First" << endl; hEff->Draw(); + hEffD->Draw("same"); // hMCPtGen->Draw(); // hMCPtPri->Draw("same"); } else { cout << "Same" << endl; hEff->Draw("same"); + hEffD->Draw("same"); // hMCPtGen->Draw(""); // hMCPtPri->Draw("same"); } @@ -408,8 +447,8 @@ void ShowAcceptanceInVzSlices() { } //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw(); - hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw(""); - hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen )->Draw("same"); + // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw(""); + // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen )->Draw("same"); } @@ -440,3 +479,21 @@ static Double_t HistoSum(const double * x, const double* p){ return value; } + +void PrintCanvas(TCanvas* c,const TString formats) { + // print a canvas in every of the given comma-separated formats + // ensure the canvas is updated + if(!doPrint) return; + c->Update(); + gSystem->ProcessEvents(); + TObjArray * arr = formats.Tokenize(","); + TIterator * iter = arr->MakeIterator(); + TObjString * element = 0; + TString name =c ->GetName(); + name.ReplaceAll(" ","_"); + name.ReplaceAll("+","Plus"); + name.ReplaceAll("-",""); + while ((element = (TObjString*) iter->Next())) { + c->Print(name+ "."+element->GetString()); + } +} diff --git a/PWG0/multPbPb/run.C b/PWG0/multPbPb/run.C index 97535daeb02..906eaea1ba7 100644 --- a/PWG0/multPbPb/run.C +++ b/PWG0/multPbPb/run.C @@ -1,12 +1,15 @@ // TODO: // 1. Check cuts for 2010 (Jochen?) // 2. Run with many centrality bins at once +#include -enum { kMyRunModeLocal = 0, kMyRunModeCAF}; +enum { kMyRunModeLocal = 0, kMyRunModeCAF, kMyRunModeGRID}; + +TList * listToLoad = new TList(); TChain * GetAnalysisChain(const char * incollection); -void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kFALSE, Int_t runMode = 0, Bool_t isMC = 0, Int_t centrBin = 0, const char * centrEstimator = "VOM", const char* option = "",Int_t workers = -1) +void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kFALSE, Int_t runMode = 0, Bool_t isMC = 0, Int_t centrBin = 0, const char * centrEstimator = "VOM", const char* option = "",TString customSuffix = "", Int_t workers = -1) { // runMode: // @@ -33,6 +36,21 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF mgr->SetMCtruthEventHandler(handler); } + + // If we are running on grid, we need the alien handler + if (runMode == kMyRunModeGRID) { + // Create and configure the alien handler plugin + gROOT->LoadMacro("CreateAlienHandler.C"); + AliAnalysisGrid *alienHandler = CreateAlienHandler(data, listToLoad, "test", isMC); + if (!alienHandler) { + cout << "Cannot create alien handler" << endl; + exit(1); + } + mgr->SetGridHandler(alienHandler); + } + + + // physics selection gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); physicsSelectionTask = AddTaskPhysicsSelection(isMC); @@ -58,7 +76,7 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF doSave = kTRUE; } - AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(kTRUE); + AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); TString pathsuffix = ""; // cuts->SetPtRange(0.15,0.2);// FIXME pt cut // const char * pathsuffix = "_pt_015_020_nofakes"; @@ -88,6 +106,11 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/multPbPb/AddTaskMultPbPbTracks.C"); AliAnalysisTaskMultPbTracks * task = AddTaskMultPbPbTracks("multPbPbtracks.root", cuts); // kTRUE enables DCA cut task->SetIsMC(useMCKinematics); + if(useMCKinematics) task->GetHistoManager()->SetSuffix("MC"); + if(customSuffix!=""){ + cout << "Setting custom suffix: " << customSuffix << endl; + task->GetHistoManager()->SetSuffix(customSuffix); + } task->SetCentralityBin(centrBin); task->SetCentralityEstimator(centrEstimator); @@ -103,10 +126,13 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF mgr->StartAnalysis("local",chain,nev); } else if (runMode == kMyRunModeCAF) { mgr->StartAnalysis("proof",TString(data)+"#esdTree",nev); + } else if (runMode == kMyRunModeGRID) { + mgr->StartAnalysis("grid"); } else { cout << "ERROR: unknown run mode" << endl; } + pathsuffix = pathsuffix + "_" + centrEstimator + "_bin_"+long(centrBin); if (doSave) MoveOutput(data, pathsuffix.Data()); @@ -159,6 +185,14 @@ TChain * GetAnalysisChain(const char * incollection){ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug=0) { + // Loads libs and par files + custom task and classes + + // Custom stuff to be loaded + listToLoad->Add(new TObjString("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx+")); + listToLoad->Add(new TObjString("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+")); + listToLoad->Add(new TObjString("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+")); + listToLoad->Add(new TObjString("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+")); + if (runMode == kMyRunModeCAF) { @@ -185,7 +219,7 @@ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug } else { - cout << "Init in Local mode" << endl; + cout << "Init in Local or Grid mode" << endl; gSystem->Load("libVMC"); gSystem->Load("libTree"); @@ -201,30 +235,16 @@ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background/")); } // Load helper classes - // TODO: replace this by a list of TOBJStrings - TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+"); - TString histoManName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+"); - TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+"); - - gSystem->ExpandPathName(taskName); - gSystem->ExpandPathName(histoManName); - gSystem->ExpandPathName(listName); - - - - // Create, add task - if (runMode == kMyRunModeCAF) { - gProof->Load(listName+(debug?"+g":"")); - gProof->Load(histoManName+(debug?"+g":"")); - gProof->Load(taskName+(debug?"+g":"")); - gProof->Load("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx++g"); - } else { - gROOT->LoadMacro(listName+(debug?"+g":"")); - gROOT->LoadMacro(histoManName+(debug?"+g":"")); - gROOT->LoadMacro(taskName+(debug?"+g":"")); - gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx++g"); - + TIterator * iter = listToLoad->MakeIterator(); + TObjString * name = 0; + while (name = (TObjString *)iter->Next()) { + gSystem->ExpandPathName(name->String()); + cout << name->String().Data(); + if (runMode == kMyRunModeCAF) { + gProof->Load(name->String()+(debug?"+g":"")); + } else { + gROOT->LoadMacro(name->String()+(debug?"+g":"")); + } } - } diff --git a/PWG0/multPbPb/run.sh b/PWG0/multPbPb/run.sh index 94d78de52ef..004fef69ae4 100755 --- a/PWG0/multPbPb/run.sh +++ b/PWG0/multPbPb/run.sh @@ -16,11 +16,12 @@ analysismode=9; #SPD + field on centrBin=0 centrEstimator="V0M" runTriggerStudy=no +customSuffix="" give_help() { cat < Number of events to be analized Misc -d Dataset or data collection (according to run mode) [default=$dataset] + - local mode: a single ESD file, an xml collection of files on + grid or a text file with a ESD per line + - caf mode: a dataset + - grid mode: a directory on alien Options specific to the multiplicity analysis -b Set centrality bin [default=$centrBin] -e Set centrality estimator [default=$centrEstimator] @@ -56,18 +62,22 @@ Available options: * == can be used in trigger studies task -t