--- /dev/null
+// ### Settings that make sense when using the Alien plugin
+//==============================================================================
+Int_t runOnData = 0; // Set to 1 if processing real data
+Int_t iCollision = 0; // 0=pp, 1=Pb-Pb
+//==============================================================================
+Bool_t usePhysicsSelection = kTRUE; // use physics selection
+Bool_t useTender = kFALSE; // use tender wagon
+Bool_t useCentrality = kTRUE; // centrality
+Bool_t useV0tender = kFALSE; // use V0 correction in tender
+Bool_t useDBG = kTRUE; // activate debugging
+Bool_t useMC = kTRUE; // use MC info
+Bool_t useKFILTER = kTRUE; // use Kinematics filter
+Bool_t useTR = kTRUE; // use track references
+Bool_t useCORRFW = kFALSE; // do not change
+Bool_t useAODTAGS = kFALSE; // use AOD tags
+Bool_t useSysInfo = kFALSE; // use sys info
+
+// ### Analysis modules to be included. Some may not be yet fully implemented.
+//==============================================================================
+Int_t iAODhandler = 1; // Analysis produces an AOD or dAOD's
+Int_t iESDMCLabelAddition= 0;
+Int_t iESDfilter = 1; // ESD to AOD filter (barrel + muon tracks)
+Int_t iMUONcopyAOD = 1; // Task that copies only muon events in a separate AOD (PWG3)
+Int_t iJETAN = 0; // Jet analysis (PWG4)
+Int_t iJETANdelta = 0; // Jet delta AODs
+Int_t iPWGHFvertexing = 0; // Vertexing HF task (PWG3)
+Int_t iPWGDQJPSIfilter = 0; // JPSI filtering (PWG3)
+Int_t iPWGHFd2h = 0; // D0->2 hadrons (PWG3)
+Bool_t doPIDResponse = 1;
+Bool_t doPIDqa = 1; //new
+
+// ### Configuration macros used for each module
+//==============================================================================
+ TString configPWGHFd2h = (iCollision==0)?"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF.C"
+ :"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF_Pb_AllCent.C";
+
+// Temporaries.
+class AliOADBPhysicsSelection;
+AliOADBPhysicsSelection *CreateOADBphysicsSelection();
+void AODmerge();
+void AddAnalysisTasks(Int_t);
+Bool_t LoadCommonLibraries();
+Bool_t LoadAnalysisLibraries();
+Bool_t LoadLibrary(const char *);
+TChain *CreateChain();
+
+//______________________________________________________________________________
+void AODtrain(Int_t merge=0)
+{
+// Main analysis train macro.
+ AliLog::SetGlobalDebugLevel(3);
+ if (merge) {
+ TGrid::Connect("alien://");
+ if (!gGrid || !gGrid->IsConnected()) {
+ ::Error("QAtrain", "No grid connection");
+ return;
+ }
+ }
+ // Set temporary merging directory to current one
+ gSystem->Setenv("TMPDIR", gSystem->pwd());
+ // Set temporary compilation directory to current one
+ gSystem->SetBuildDir(gSystem->pwd(), kTRUE);
+ printf("==================================================================\n");
+ printf("=========== RUNNING FILTERING TRAIN ==========\n");
+ printf("==================================================================\n");
+ printf("= Configuring analysis train for: =\n");
+ if (usePhysicsSelection) printf("= Physics selection =\n");
+ if (useTender) printf("= TENDER =\n");
+ if (iESDfilter) printf("= ESD filter =\n");
+ if (iMUONcopyAOD) printf("= MUON copy AOD =\n");
+ if (iJETAN) printf("= Jet analysis =\n");
+ if (iJETANdelta) printf("= Jet delta AODs =\n");
+ if (iPWGHFvertexing) printf("= PWGHF vertexing =\n");
+ if (iPWGDQJPSIfilter) printf("= PWGDQ j/psi filter =\n");
+ if (iPWGHFd2h) printf("= PWGHF D0->2 hadrons QA =\n");
+
+ // Load common libraries and set include path
+ if (!LoadCommonLibraries()) {
+ ::Error("AnalysisTrain", "Could not load common libraries");
+ return;
+ }
+
+ // Make the analysis manager and connect event handlers
+ AliAnalysisManager *mgr = new AliAnalysisManager("Analysis Train", "Production train");
+ if (useSysInfo) mgr->SetNSysInfo(100);
+ // Load analysis specific libraries
+ if (!LoadAnalysisLibraries()) {
+ ::Error("AnalysisTrain", "Could not load analysis libraries");
+ return;
+ }
+
+ // Create input handler (input container created automatically)
+ // ESD input handler
+ AliESDInputHandler *esdHandler = new AliESDInputHandler();
+ mgr->SetInputEventHandler(esdHandler);
+ // Monte Carlo handler
+ if (useMC) {
+ AliMCEventHandler* mcHandler = new AliMCEventHandler();
+ mgr->SetMCtruthEventHandler(mcHandler);
+ mcHandler->SetReadTR(useTR);
+ }
+ // AOD output container, created automatically when setting an AOD handler
+ if (iAODhandler) {
+ // AOD output handler
+ AliAODHandler* aodHandler = new AliAODHandler();
+ aodHandler->SetOutputFileName("AliAOD.root");
+ mgr->SetOutputEventHandler(aodHandler);
+ }
+ // Debugging if needed
+ if (useDBG) mgr->SetDebugLevel(3);
+
+ AddAnalysisTasks(merge);
+ if (merge) {
+ AODmerge();
+ mgr->InitAnalysis();
+ mgr->SetGridHandler(new AliAnalysisAlien);
+ mgr->StartAnalysis("gridterminate",0);
+ return;
+ }
+ // Run the analysis
+ //
+ TChain *chain = CreateChain();
+ if (!chain) return;
+
+ TStopwatch timer;
+ timer.Start();
+ mgr->SetSkipTerminate(kTRUE);
+ if (mgr->InitAnalysis()) {
+ mgr->PrintStatus();
+ mgr->StartAnalysis("local", chain);
+ }
+ timer.Print();
+}
+
+//______________________________________________________________________________
+void AddAnalysisTasks(Int_t merge){
+ // Add all analysis task wagons to the train
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+
+ //
+ // Tender and supplies. Needs to be called for every event.
+ //
+ AliAnalysisManager::SetCommonFileName("AODQA.root");
+ if (useTender) {
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/TenderSupplies/AddTaskTender.C");
+ // IF V0 tender needed, put kTRUE below
+ AliAnalysisTaskSE *tender = AddTaskTender(useV0tender);
+// tender->SetDebugLevel(2);
+ }
+
+ if (usePhysicsSelection) {
+ // Physics selection task
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ mgr->RegisterExtraFile("event_stat.root");
+ AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(kFALSE);
+// AliOADBPhysicsSelection * oadbDefaultPbPb = CreateOADBphysicsSelection();
+// physSelTask->GetPhysicsSelection()->SetCustomOADBObjects(oadbDefaultPbPb,0,0);
+// if (!merge) mgr->AddStatisticsTask(AliVEvent::kAny);
+ }
+ // Centrality (only Pb-Pb)
+ if (iCollision && useCentrality) {
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+ AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
+ taskCentrality->SelectCollisionCandidates(AliVEvent::kAny);
+ }
+
+ if(iESDMCLabelAddition) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWG/muondep/AddTaskESDMCLabelAddition.C");
+ AliAnalysisTaskESDMCLabelAddition *esdmclabel = AddTaskESDMCLabelAddition(kFALSE);
+ }
+
+ if (iESDfilter) {
+ // ESD filter task configuration.
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
+ if (iMUONcopyAOD) {
+ printf("Registering delta AOD file\n");
+ mgr->RegisterExtraFile("AliAOD.Muons.root");
+ mgr->RegisterExtraFile("AliAOD.Dimuons.root");
+ AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kTRUE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,kTRUE,kTRUE,kTRUE,1100,1);
+ } else {
+ AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kFALSE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,kTRUE,kTRUE,kTRUE,1100,1);
+ }
+ }
+
+// ********** PWG3 wagons ******************************************************
+ // PWGHF vertexing
+ if (iPWGHFvertexing) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/AddTaskVertexingHF.C");
+ if (!iPWGHFd2h) TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
+ AliAnalysisTaskSEVertexingHF *taskvertexingHF = AddTaskVertexingHF();
+ if (!taskvertexingHF) ::Warning("AnalysisTrainNew", "AliAnalysisTaskSEVertexingHF cannot run for this train conditions - EXCLUDED");
+ else mgr->RegisterExtraFile("AliAOD.VertexingHF.root");
+ taskvertexingHF->SelectCollisionCandidates(0);
+ }
+
+ // PWGDQ JPSI filtering (only pp)
+ if (iPWGDQJPSIfilter && (iCollision==0)) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWGDQ/dielectron/macros/AddTaskJPSIFilter.C");
+ AliAnalysisTaskSE *taskJPSIfilter = AddTaskJPSIFilter();
+ if (!taskJPSIfilter) ::Warning("AnalysisTrainNew", "AliAnalysisTaskDielectronFilter cannot run for this train conditions - EXCLUDED");
+ else mgr->RegisterExtraFile("AliAOD.Dielectron.root");
+ taskJPSIfilter->SelectCollisionCandidates(0);
+ }
+
+ // PWGHF D2h
+ if (iPWGHFd2h) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/AddD2HTrain.C");
+ TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
+ AddD2HTrain(kFALSE, 1,0,0,0,0,0,0,0,0,0,0);
+ }
+
+ // ********** PWG4 wagons ******************************************************
+ // Jet analysis
+
+ // Configurations flags, move up?
+ TString kDeltaAODJetName = "AliAOD.Jets.root"; //
+ Bool_t kIsPbPb = (iCollision==0)?false:true; // can be more intlligent checking the name of the data set
+ TString kDefaultJetBackgroundBranch = "";
+ TString kJetSubtractBranches = "";
+ UInt_t kHighPtFilterMask = 128;// from esd filter
+ UInt_t iPhysicsSelectionFlag = AliVEvent::kMB;
+ if (iJETAN) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJets.C");
+ // Default jet reconstructor running on ESD's
+ AliAnalysisTaskJets *taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,0); // no background subtraction
+ if (!taskjets) ::Fatal("AnalysisTrainNew", "AliAnalysisTaskJets cannot run for this train conditions - EXCLUDED");
+ if(kDeltaAODJetName.Length()>0) taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
+ if (iJETANdelta) {
+ // AddTaskJetsDelta("AliAOD.Jets.root"); // need to modify this accordingly in the add task jets
+ mgr->RegisterExtraFile(kDeltaAODJetName.Data());
+ TString cTmp("");
+ if(kIsPbPb){
+ // UA1 intrinsic background subtraction
+ taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,2); // background subtraction
+ if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
+ }
+ // SICONE
+ taskjets = AddTaskJets("AOD","SISCONE",0.4,kHighPtFilterMask,0.15,0); //no background subtration to be done later....
+ if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
+ cTmp = taskjets->GetNonStdBranch();
+ if(cTmp.Length()>0)kJetSubtractBranches += Form("%s ",cTmp.Data());
+
+ // Add the clusters..
+ gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJetCluster.C");
+ AliAnalysisTaskJetCluster *taskCl = 0;
+ Float_t fCenUp = 0;
+ Float_t fCenLo = 0;
+ Float_t fTrackEtaWindow = 0.9;
+ taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"KT",0.4,0,1, kDeltaAODJetName.Data(),0.15,fTrackEtaWindow,0); // this one is for the background and random jets, random cones with no skip
+ taskCl->SetBackgroundCalc(kTRUE);
+ taskCl->SetNRandomCones(10);
+ taskCl->SetCentralityCut(fCenLo,fCenUp);
+ taskCl->SetGhostEtamax(fTrackEtaWindow);
+ kDefaultJetBackgroundBranch = Form("%s_%s",AliAODJetEventBackground::StdBranchName(),taskCl->GetJetOutputBranch());
+
+ taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.4,2,1,kDeltaAODJetName.Data(),0.15);
+ taskCl->SetCentralityCut(fCenLo,fCenUp);
+ if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
+ taskCl->SetNRandomCones(10);
+ kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
+
+ taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.2,0,1,kDeltaAODJetName.Data(),0.15);
+ taskCl->SetCentralityCut(fCenLo,fCenUp);
+ if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
+ kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
+
+ // DO THE BACKGROUND SUBTRACTION
+ if(kIsPbPb&&kJetSubtractBranches.Length()){
+ gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJetBackgroundSubtract.C");
+ AliAnalysisTaskJetBackgroundSubtract *taskSubtract = 0;
+ taskSubtract = AddTaskJetBackgroundSubtract(kJetSubtractBranches,1,"B0","B%d");
+ taskSubtract->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
+ if(kDeltaAODJetName.Length()>0)taskSubtract->SetNonStdOutputFile(kDeltaAODJetName.Data());
+ }
+ }
+ }
+}
+
+//______________________________________________________________________________
+Bool_t LoadCommonLibraries()
+{
+// Load common analysis libraries.
+ if (!gSystem->Getenv("ALICE_ROOT")) {
+ ::Error("AnalysisTrainNew.C::LoadCommonLibraries", "Analysis train requires that analysis libraries are compiled with a local AliRoot");
+ return kFALSE;
+ }
+ Bool_t success = kTRUE;
+ // Load framework classes. Par option ignored here.
+ success &= LoadLibrary("libSTEERBase.so");
+ success &= LoadLibrary("libESD.so");
+ success &= LoadLibrary("libAOD.so");
+ success &= LoadLibrary("libANALYSIS.so");
+ success &= LoadLibrary("libOADB.so");
+ success &= LoadLibrary("libANALYSISalice.so");
+ success &= LoadLibrary("libCORRFW.so");
+ gROOT->ProcessLine(".include $ALICE_ROOT/include");
+ if (success) {
+ ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries: SUCCESS");
+ ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Include path for Aclic compilation:\n%s",
+ gSystem->GetIncludePath());
+ } else {
+ ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries: FAILED");
+ }
+ return success;
+}
+
+//______________________________________________________________________________
+Bool_t LoadAnalysisLibraries()
+{
+// Load common analysis libraries.
+ if (useTender) {
+ if (!LoadLibrary("TENDER") ||
+ !LoadLibrary("TENDERSupplies")) return kFALSE;
+ }
+ if (iESDfilter || iPWGMuonTrain) {
+ if (!LoadLibrary("PWGmuon")) return kFALSE;
+ }
+ if (iESDMCLabelAddition) {
+ if (!LoadLibrary("PWGmuondep")) return kFALSE;
+ }
+ // JETAN
+ if (iJETAN) {
+ if (!LoadLibrary("JETAN")) return kFALSE;
+ }
+ if (iJETANdelta) {
+ if (!LoadLibrary("JETAN") ||
+ !LoadLibrary("CGAL") ||
+ !LoadLibrary("fastjet") ||
+ !LoadLibrary("siscone") ||
+ !LoadLibrary("SISConePlugin") ||
+ !LoadLibrary("FASTJETAN")) return kFALSE;
+ }
+ // PWG3 Vertexing HF
+ if (iPWGHFvertexing || iPWGHFd2h) {
+ if (!LoadLibrary("PWGflowBase") ||
+ !LoadLibrary("PWGflowTasks") ||
+ !LoadLibrary("PWGHFvertexingHF")) return kFALSE;
+ }
+ // if (iPWGHFvertexing || iPWG3d2h) {
+ // if (!LoadLibrary("PWG3base") ||
+ // !LoadLibrary("PWGHFvertexingHF")) return kFALSE;
+ // }
+ // PWG3 dielectron
+ if (iPWGDQJPSIfilter) {
+ if (!LoadLibrary("PWGDQdielectron")) return kFALSE;
+ }
+
+ ::Info("AnalysisTrainNew.C::LoadAnalysisLibraries", "Load other libraries: SUCCESS");
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+Bool_t LoadLibrary(const char *module)
+{
+// Load a module library in a given mode. Reports success.
+ Int_t result;
+ TString mod(module);
+ if (!mod.Length()) {
+ ::Error("AnalysisTrainNew.C::LoadLibrary", "Empty module name");
+ return kFALSE;
+ }
+ // If a library is specified, just load it
+ if (mod.EndsWith(".so")) {
+ mod.Remove(mod.Index(".so"));
+ result = gSystem->Load(mod);
+ if (result < 0) {
+ ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load library %s", module);
+ return kFALSE;
+ }
+ return kTRUE;
+ }
+ // Check if the library is already loaded
+ if (strlen(gSystem->GetLibraries(Form("%s.so", module), "", kFALSE)) > 0) return kTRUE;
+ result = gSystem->Load(Form("lib%s.so", module));
+ if (result < 0) {
+ ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load module %s", module);
+ return kFALSE;
+ }
+ return kTRUE;
+}
+
+
+//______________________________________________________________________________
+TChain *CreateChain()
+{
+// Create the input chain
+ chain = new TChain("esdTree");
+ if (gSystem->AccessPathName("AliESDs.root"))
+ ::Error("AnalysisTrainNew.C::CreateChain", "File: AliESDs.root not in ./data dir");
+ else
+ chain->Add("AliESDs.root");
+ if (chain->GetNtrees()) return chain;
+ return NULL;
+}
+
+//______________________________________________________________________________
+void AODmerge()
+{
+// Merging method. No staging and no terminate phase.
+ TStopwatch timer;
+ timer.Start();
+ TString outputDir = "wn.xml";
+ TString outputFiles = "EventStat_temp.root,AliAOD.root,AliAOD.Muons.root";
+ TString mergeExcludes = "";
+ TObjArray *list = outputFiles.Tokenize(",");
+ TIter *iter = new TIter(list);
+ TObjString *str;
+ TString outputFile;
+ Bool_t merged = kTRUE;
+ while((str=(TObjString*)iter->Next())) {
+ outputFile = str->GetString();
+ // Skip already merged outputs
+ if (!gSystem->AccessPathName(outputFile)) {
+ printf("Output file <%s> found. Not merging again.",outputFile.Data());
+ continue;
+ }
+ if (mergeExcludes.Contains(outputFile.Data())) continue;
+ merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, 0);
+ if (!merged) {
+ printf("ERROR: Cannot merge %s\n", outputFile.Data());
+ return;
+ }
+ }
+ // all outputs merged, validate
+ ofstream out;
+ out.open("outputs_valid", ios::out);
+ out.close();
+ timer.Print();
+}
fCurrentDigit(0),
fCurrentCluster(0),
fSegmentation(0),
- fNPlanes(0)
+ fNPlanes(0),
+ sw(0)
{
// Default constructor
for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) fClustersPerPlane[iPlane] = NULL;
fDigitsInCluster = new TClonesArray("AliMFTDigit", fNMaxDigitsPerCluster);
fDigitsInCluster -> SetOwner(kTRUE);
+ sw = new TStopwatch();
+ sw -> Reset();
}
}
if (fDigitsInCluster) fDigitsInCluster->Delete(); delete fDigitsInCluster; fDigitsInCluster = NULL;
delete fSegmentation; fSegmentation=NULL;
+
+ delete sw;
AliDebug(1, "... done!");
AliInfo("Starting Clusterization for MFT");
AliDebug(1, Form("nPlanes = %d", fNPlanes));
+ if (!sw) sw = new TStopwatch();
+ sw -> Reset();
+
StartEvent();
Bool_t isDigAvailableForNewCluster = kTRUE;
AliDebug(1, Form("Plane %02d", iPlane));
+ Int_t nDetElem = fSegmentation->GetPlane(iPlane)->GetNActiveElements();
+ TClonesArray *clustersPerDetElem[fNMaxDetElemPerPlane] = {0};
+ for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) clustersPerDetElem[iDetElem] = new TClonesArray("AliMFTCluster");
+
myDigitList = (TClonesArray*) pDigitList->At(iPlane);
AliDebug(1, Form("myDigitList->GetEntries() = %d", myDigitList->GetEntries()));
Int_t cycleOverDigits = 0;
- Double_t myCutForAvailableDigits = fCutForAvailableDigits;
+ Double_t myCutForAvailableDigits = 0;
+ Int_t currentDetElem = -1;
+ Int_t currentDetElemLocal = -1;
+ Bool_t areThereSkippedDigits = kFALSE;
+
+ sw -> Start();
+
while (myDigitList->GetEntries()) {
for (Int_t iDig=0; iDig<myDigitList->GetEntries(); iDig++) {
- AliDebug(1, Form("Check %d: Digit %5d of %5d", cycleOverDigits, iDig, myDigitList->GetEntries()));
-
fCurrentDigit = (AliMFTDigit*) myDigitList->At(iDig);
+
+ if (!iDig) {
+ if (fCurrentDigit->GetDetElemID() != currentDetElem) {
+ // first iteration over the digits of a specific detection element
+ currentDetElem = fCurrentDigit->GetDetElemID();
+ currentDetElemLocal = fSegmentation->GetDetElemLocalID(currentDetElem);
+ // printf("Reading a new detection element: %d\n",currentDetElem);
+ cycleOverDigits = 0;
+ myCutForAvailableDigits = fCutForAvailableDigits;
+ }
+ else if (fCurrentDigit->GetDetElemID()==currentDetElem && areThereSkippedDigits) {
+ // second (or further) iteration over the digits of a specific detection element
+ cycleOverDigits++;
+ myCutForAvailableDigits -= 0.5;
+ }
+ areThereSkippedDigits = kFALSE;
+ }
+ else {
+ areThereSkippedDigits = kTRUE;
+ if (fCurrentDigit->GetDetElemID() != currentDetElem) break;
+ }
+
isDigAvailableForNewCluster = kTRUE;
- for (Int_t iCluster=0; iCluster<fClustersPerPlane[iPlane]->GetEntries(); iCluster++) {
- fCurrentCluster = (AliMFTCluster*) fClustersPerPlane[iPlane]->At(iCluster);
- AliDebug(2, Form("Distance between cluster and digit = %f",fCurrentCluster->GetDistanceFromPixel(fCurrentDigit)));
+ // AliDebug(2, Form("Check %d: Digit %5d of %5d", cycleOverDigits, iDig, myDigitList->GetEntries()));
+
+ for (Int_t iCluster=0; iCluster<clustersPerDetElem[currentDetElemLocal]->GetEntries(); iCluster++) {
+ fCurrentCluster = (AliMFTCluster*) clustersPerDetElem[currentDetElemLocal]->At(iCluster);
+ // AliDebug(2, Form("Distance between cluster and digit = %f",fCurrentCluster->GetDistanceFromPixel(fCurrentDigit)));
if (fCurrentCluster->GetDistanceFromPixel(fCurrentDigit) < fCutForAttachingDigits) {
fCurrentCluster->AddPixel(fCurrentDigit);
myDigitList->Remove(fCurrentDigit);
myDigitList->Remove(fCurrentDigit);
myDigitList->Compress();
iDig--;
- new ((*fClustersPerPlane[iPlane])[fClustersPerPlane[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
- delete newCluster;
+ new ((*clustersPerDetElem[currentDetElemLocal])[clustersPerDetElem[currentDetElemLocal]->GetEntries()]) AliMFTCluster(*newCluster);
+ delete newCluster;
}
-
+
} // end of cycle over the digits
- if (cycleOverDigits) myCutForAvailableDigits -= 0.5;
- cycleOverDigits++;
-
} // no more digits to check in current plane!
- for (Int_t iCluster=0; iCluster<fClustersPerPlane[iPlane]->GetEntries(); iCluster++) {
- fCurrentCluster = (AliMFTCluster*) fClustersPerPlane[iPlane]->At(iCluster);
- fCurrentCluster -> TerminateCluster();
+ sw -> Print("m");
+
+ printf("Plane %d: clusters found in %f seconds\n",iPlane,sw->CpuTime());
+ sw->Start();
+
+ // Now we merge the cluster lists coming from each detection element, to build the cluster list of the plane
+
+ AliMFTCluster *newCluster = NULL;
+ for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) {
+ for (Int_t iCluster=0; iCluster<clustersPerDetElem[iDetElem]->GetEntries(); iCluster++) {
+ newCluster = (AliMFTCluster*) (clustersPerDetElem[iDetElem]->At(iCluster));
+ newCluster -> TerminateCluster();
+ // new ((*fClustersPerPlane[iPlane])[fClustersPerPlane[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
+ }
}
- AliDebug(1, Form("Found %d clusters in plane %02d", fClustersPerPlane[iPlane]->GetEntries(), iPlane));
+ printf("%d Clusters in plane %02d merged in %f seconds\n", fClustersPerPlane[iPlane]->GetEntries(), iPlane, sw->CpuTime());
+
+// sw->Start();
+// Double_t centerX=0, centerY=0;
+// for (Int_t iCluster=0; iCluster<fClustersPerPlane[iPlane]->GetEntries(); iCluster++) {
+// newCluster = (AliMFTCluster*) (fClustersPerPlane[iPlane]->At(iCluster));
+// centerX += newCluster-> GetX();
+// centerY += newCluster-> GetY();
+// }
+// centerX /= Double_t(fClustersPerPlane[iPlane]->GetEntries());
+// centerY /= Double_t(fClustersPerPlane[iPlane]->GetEntries());
+// Double_t time = sw->CpuTime();
+// printf("Loop over %d performed in %f seconds. <X> = %f, <Y> = %f\n",fClustersPerPlane[iPlane]->GetEntries(),time,centerX,centerY);
+
+
+ for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) {
+ clustersPerDetElem[iDetElem] -> Delete();
+ delete clustersPerDetElem[iDetElem];
+ }
myDigitList -> Delete();
fSegmentation(0),
fNPlanesMFT(0),
fNPlanesMFTAnalyzed(0),
- fSigmaClusterCut(0),
+ fSigmaClusterCut(2),
fScaleSigmaClusterCut(1.),
fNMaxMissingMFTClusters(0),
fGlobalTrackingDiverged(kFALSE),
SetNPlanesMFT(fSegmentation->GetNPlanes());
AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations
+ AliInfo(Form("fMFT = %p, fSegmentation = %p", fMFT, fSegmentation));
+
for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
- fMFTClusterArray[iPlane] = 0;
+ fMFTClusterArray[iPlane] = new TClonesArray("AliMFTCluster");
fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
fMFTClusterArrayBack[iPlane] = new TClonesArray("AliMFTCluster");
- fIsPlaneMandatory[iPlane] = kFALSE;
+ fMFTClusterArray[iPlane] -> SetOwner(kTRUE);
+ fMFTClusterArrayFront[iPlane] -> SetOwner(kTRUE);
+ fMFTClusterArrayBack[iPlane] -> SetOwner(kTRUE);
+ fMinResearchRadiusAtPlane[iPlane] = 0.;
}
+ fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
+
}
//====================================================================================================================================================
// destructor
for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+ fMFTClusterArray[iPlane] -> Delete();
delete fMFTClusterArray[iPlane];
delete fMFTClusterArrayFront[iPlane];
delete fMFTClusterArrayBack[iPlane];
}
+ delete fCandidateTracks;
+
}
//====================================================================================================================================================
// This function loads the MFT clusters
//--------------------------------------------------------------------
+ for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+ AliDebug(1, Form("Setting Address for Branch Plane_%02d", iPlane));
+ cTree->SetBranchAddress(Form("Plane_%02d",iPlane), &fMFTClusterArray[iPlane]);
+ }
+
if (!cTree->GetEvent()) return kFALSE;
for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
- AliDebug(1, Form("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries()));
- fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
+ AliInfo(Form("plane %02d: nClusters = %d\n", iPlane, fMFTClusterArray[iPlane]->GetEntries()));
}
SeparateFrontBackClusters();
void AliMFTTrackerMU::UnloadClusters() {
- //--------------------------------------------------------------------
- // This function unloads MFT clusters
- //--------------------------------------------------------------------
+// //--------------------------------------------------------------------
+// // This function unloads MFT clusters
+// //--------------------------------------------------------------------
for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
fMFTClusterArray[iPlane] -> Clear("C");
// The clusters must be already loaded !
//--------------------------------------------------------------------
+ // Tree of AliMuonForwardTrack objects. Created outside the ESD framework for cross-check purposes
+
+ TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
TTree *outputTreeMuonGlobalTracks = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
TClonesArray *muonForwardTracks = new TClonesArray("AliMuonForwardTrack");
outputTreeMuonGlobalTracks -> Branch("tracks", &muonForwardTracks);
fESD = event;
+ // get the ITS primary vertex
+
+ Double_t vertex[3] = {0., 0., 0.};
+ const AliESDVertex* esdVert = fESD->GetVertex();
+ if (esdVert->GetNContributors() > 0 || !strcmp(esdVert->GetTitle(),"vertexer: smearMC")) {
+ esdVert->GetXYZ(vertex);
+ AliDebug(1,Form("found vertex (%e,%e,%e)",vertex[0],vertex[1],vertex[2]));
+ }
+
//----------- Read ESD MUON tracks -------------------
Int_t nTracksMUON = event->GetNumberOfMuonTracks();
fNPlanesMFTAnalyzed = 0;
- AliDebug(1, "**************************************************************************************\n");
- AliDebug(1, Form("*************************** MUON TRACK %3d ***************************************\n", iTrack));
- AliDebug(1, "**************************************************************************************\n");
+ AliInfo("****************************************************************************************");
+ AliInfo(Form("*************************** MUON TRACK %3d/%d ***************************************", iTrack, nTracksMUON));
+ AliInfo("****************************************************************************************");
fCandidateTracks -> Delete();
fNPlanesMFTAnalyzed = 0;
const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
- fMUONTrack = NULL;
- AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
+ if (fMUONTrack) delete fMUONTrack;
+ fMUONTrack = new AliMUONTrack();
+
+ AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kTRUE); // last arguments should be kTRUE or kFALSE??
// the track we are going to build, starting from fMUONTrack and adding the MFT clusters
AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
Int_t nCandidates = fCandidateTracks->GetEntriesFast();
for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iCandidate);
+
// if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
// (several new tracks can be created for one old track)
if (FindClusterInPlane(iPlane) == kDiverged) {
fCandidateTracks->Compress();
}
-
+
// -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
fGlobalTrackingDiverged = kFALSE;
// If we have several final tracks, we must find the best candidate:
Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
- AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
-
+ AliInfo(Form("nFinalTracks = %d", nFinalTracks));
+
Int_t nGoodClustersBestCandidate = 0;
Int_t idBestCandidate = 0;
Double_t bestChi2 = -1.; // variable defining the best candidate
AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
- //----------------------- Save the information to the AliESDMuonForwardTrack object
+ new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
-// AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
-// myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
-// myESDTrack -> SetChi2(newTrack->GetGlobalChi2());
-// myESDTrack -> SetCharge(newTrack->GetCharge());
-// myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
-
- //---------------------------------------------------------------------------------
+ //----------------------- Save the information to the AliESDMuonGlobalTrack object
+
+ newTrack -> EvalKinem(vertex[2]); // we evaluate the kinematics at the primary vertex
+
+ AliDebug(2,"Creating a new Muon Global Track");
+ AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
+ myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
+ myESDTrack -> SetChi2OverNdf(newTrack->GetChi2OverNdf());
+ myESDTrack -> SetCharge(newTrack->GetCharge());
+ myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
+ myESDTrack -> SetFirstTrackingPoint(newTrack->GetMFTCluster(0)->GetX(), newTrack->GetMFTCluster(0)->GetY(), newTrack->GetMFTCluster(0)->GetZ());
+ myESDTrack -> SetXYAtVertex(newTrack->GetOffsetX(vertex[0], vertex[2]), newTrack->GetOffsetX(vertex[1], vertex[2]));
+ myESDTrack -> SetRAtAbsorberEnd(newTrack->GetRAtAbsorberEnd());
+ myESDTrack -> SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
+ myESDTrack -> SetMuonClusterMap(esdTrack->GetMuonClusterMap());
+ myESDTrack -> SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
+ myESDTrack -> SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
+ myESDTrack -> Connected(esdTrack->IsConnected());
- new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
+ //---------------------------------------------------------------------------------
}
- fCandidateTracks->Delete();
fFinalBestCandidate = NULL;
}
+ outputTreeMuonGlobalTracks -> Fill();
+
Int_t myEventID = 0;
- TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
while (outputFileMuonGlobalTracks->cd(Form("Event%d",myEventID))) myEventID++;
outputFileMuonGlobalTracks -> mkdir(Form("Event%d",myEventID));
outputFileMuonGlobalTracks -> cd(Form("Event%d",myEventID));
outputTreeMuonGlobalTracks -> Write();
outputFileMuonGlobalTracks -> Close();
+ muonForwardTracks -> Delete();
+ delete muonForwardTracks;
+
return 0;
}
Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) {
- AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
-
// !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
// propagate track to plane #planeId (both to front and back active sensors)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff