]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/macros/electrons/anaJete.C
protect against invoking MC Handler for AOD
[u/mrichter/AliRoot.git] / PWG4 / macros / electrons / anaJete.C
CommitLineData
b21a4af6 1/* $Id: $ */\r
2//--------------------------------------------------\r
3// Example macro to do multi-platform JETAN/FASTJET analysis\r
4// Can be executed with Root and AliRoot\r
5//\r
6// Configured by options and definitions set in the lines below and\r
7// additional external configuration files and environement variables.\r
8//\r
9// Author: K. Read\r
10//\r
11//-------------------------------------------------\r
12enum anaModes {mLocal, mLocalCAF, mPROOF, mGRID, mPLUGIN};\r
13//mLocal = 0: Analyze locally files in your computer\r
14//mLocalCAF = 1: Analyze locally CAF files\r
15//mPROOF = 2: Analyze CAF files with PROOF\r
16//mGRID = 3: Analyze files on GRID\r
17//mPLUGIN = 4: Analyze files on GRID with AliEn plugin\r
18\r
19//---------------------------------------------------------------------------\r
20//Settings to read locally several files, only for "mLocal" mode\r
21//The different values are default, they can be set with environmental \r
22//variables: INDIR, PATTERN, NFILES, respectively\r
23//char * kInDir = "/afs/cern.ch/user/k/kread/public/data/"; \r
24char * kInDir = "/user/data/files";\r
25char * kPattern = ""; // Data are in files kInDir/kPattern+i\r
26Int_t kFile = 1; // Number of files\r
27//---------------------------------------------------------------------------\r
28//Collection file for grid analysis\r
29char * kXML = "collection.xml";\r
30//---------------------------------------------------------------------------\r
31//Data directory for PROOF analysis\r
32char * kmydataset = "/COMMON/COMMON/LHC09a4_run8101X";\r
33//char * kmydataset = "/PWG4/mcosenti/LHC08d10_ppElectronB_Jets#esdTree";\r
34//---------------------------------------------------------------------------\r
35//Scale histograms from file. Change to kTRUE when xsection file exists\r
36//Put name of file containing xsection \r
37//Put number of events per ESD file\r
38//This is an specific case for normalization of Pythia files.\r
9da13925 39const Bool_t kGetXSectionFromFileAndScale = kTRUE;\r
b21a4af6 40const char * kXSFileName = "pyxsec.root";\r
41const Int_t kNumberOfEventsPerFile = 200; \r
42//---------------------------------------------------------------------------\r
43\r
44const Bool_t kMC = kTRUE; //With real data kMC = kFALSE\r
d4df5eeb 45TString kInputData = "ESD";//ESD, AOD, MC\r
b21a4af6 46TString kTreeName = "esdTree";\r
47//const Bool_t kMergeAODs = kTRUE; //uncomment for AOD merging\r
48const Bool_t kMergeAODs = kFALSE; //uncomment for no AOD merging\r
d4df5eeb 49const Bool_t kUsePAR = kTRUE; //set to kFALSE for libraries\r
50const Bool_t kDoJetTask = kTRUE; //set to kFALSE to skip JETAN task\r
b21a4af6 51Int_t sevent = 0;\r
d4df5eeb 52\r
b21a4af6 53Int_t mode = mLocal;\r
54char sconfig1[1024] = "ConfigPWG4AODtoAOD"; //"ConfigAnalysis";\r
55char sconfig2[1024] = "ConfigJetAnalysisFastJet.C";//"ConfigAnalysis";\r
56char sconfig3[1024] = "ConfigAnalysisElectron"; //"ConfigAnalysis";\r
57\r
ef2c838f 58//Initialize the cross section and ntrials values. Do not modify.\r
59Double_t xsection = 0;\r
60Float_t ntrials = 0;\r
61\r
b21a4af6 62void anaJete()\r
63{\r
64 // Main\r
65\r
66 //Process environmental variables from command line:\r
67 ProcessEnvironment(); \r
d4df5eeb 68 printf("Final Variables: kInputData %s, mode %d, config2 %s, config3 %s, sevent %d\n",kInputData.Data(),mode,sconfig2,sconfig3,sevent);\r
b21a4af6 69\r
70 //--------------------------------------------------------------------\r
71 // Load analysis libraries\r
72 // Look at the method below, \r
73 // change whatever you need for your analysis case\r
74 // ------------------------------------------------------------------\r
75 LoadLibraries(mode) ;\r
76 \r
77 //-------------------------------------------------------------------------------------------------\r
78 //Create chain from ESD and from cross sections files, look below for options.\r
79 //------------------------------------------------------------------------------------------------- \r
80 if(kInputData == "ESD") kTreeName = "esdTree" ;\r
81 else if(kInputData == "AOD") kTreeName = "aodTree" ;\r
82 else if (kInputData == "MC") kTreeName = "TE" ;\r
83 else {\r
84 cout<<"Wrong data type "<<kInputData<<endl;\r
85 break;\r
86 }\r
87\r
88 TChain * chain = new TChain(kTreeName) ;\r
89 TChain * chainxs = new TChain("Xsection") ;\r
90\r
91 if (mode==mLocal || mode==mLocalCAF || mode == mGRID) {\r
92 CreateChain(mode, chain, chainxs); \r
ef2c838f 93 cout<<"Chain created"<<endl;\r
b21a4af6 94 }\r
95\r
96 if( chain || mode==mPROOF || mode==mPLUGIN ){\r
97 AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen\r
98 \r
99 //--------------------------------------\r
100 // Make the analysis manager\r
101 //-------------------------------------\r
102 AliAnalysisManager *mgr = new AliAnalysisManager("Jet Manager", "Jet Manager");\r
103\r
104 if( mode == mPLUGIN ){\r
105 // Create and configure the alien handler plugin\r
106 if (!AliAnalysisGrid::CreateToken()) return NULL;\r
107 AliAnalysisAlien *plugin = new AliAnalysisAlien();\r
108 plugin->SetRunMode("submit");\r
109 //Uncomment the following 3 lines to permit auto xml creation\r
110 //plugin->SetGridDataDir("/alice/sim/PDC_08b/LHC08d10/"); //dummy\r
111 //plugin->SetDataPattern("AliESDs.root"); //dummy\r
112 //plugin->AddRunNumber(30010); //dummy\r
113 plugin->AddDataFile("mycollect.xml");\r
114 plugin->SetGridWorkingDir("work3");\r
115 plugin->SetAdditionalLibs("anaJet.C ConfigJetAnalysisFastJet.C ConfigAnalysisElectron.C ANALYSIS.par ANALYSISalice.par AOD.par ESD.par STEERBase.par JETAN.par FASTJETAN.par");\r
116 plugin->SetJDLName("anaJet.jdl");\r
117 plugin->SetExecutable("anaJet.sh");\r
118 plugin->SetOutputFiles("histos.root");\r
119 AliAnalysisGrid *alienHandler = plugin;\r
120 if (!alienHandler) return;\r
121\r
122 // Connect plug-in to the analysis manager\r
123 mgr->SetGridHandler(alienHandler);\r
124 }\r
125\r
126 // MC handler\r
9da13925 127 if( (kMC && (kInputData == "ESD")) || kInputData == "MC"){\r
b21a4af6 128 AliMCEventHandler* mcHandler = new AliMCEventHandler();\r
129 mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file\r
130 mgr->SetMCtruthEventHandler(mcHandler);\r
131 if( kInputData == "MC") mgr->SetInputEventHandler(NULL);\r
132 }\r
133\r
134 // AOD output handler\r
135 AliAODHandler* aodoutHandler = new AliAODHandler();\r
136 aodoutHandler->SetOutputFileName("aodoutput.root");\r
137 if(kMergeAODs)aodoutHandler->SetCreateNonStandardAOD();\r
138 mgr->SetOutputEventHandler(aodoutHandler);\r
139\r
140 //input\r
141 if(kInputData == "ESD"){\r
142 // ESD handler\r
143 AliESDInputHandler *esdHandler = new AliESDInputHandler();\r
144 mgr->SetInputEventHandler(esdHandler);\r
145 }\r
146 if(kInputData == "AOD"){\r
147 // AOD handler\r
148 AliAODInputHandler *aodHandler = new AliAODInputHandler();\r
149 mgr->SetInputEventHandler(aodHandler);\r
150 if(kMergeAODs){\r
151 char path[1024];\r
152 sprintf(path,"AliAOD.root");\r
153 if(gSystem->Getenv("SIMPATH"))\r
154 sprintf(path,"%s/AliAOD.root",gSystem->Getenv("SIMPATH"));\r
155 cout<<"Config: Second input file: "<<path<<endl;\r
156 aodHandler->SetMergeEvents(kTRUE);\r
157 aodHandler->AddFriend(path);\r
158 cout<<"Config: Starting event for second input file: "<<sevent<<endl;\r
159 aodHandler->SetMergeOffset(sevent);\r
160 }\r
161 }\r
162\r
163 mgr->SetDebugLevel(3); // For debugging, do not uncomment if you want no messages.\r
164\r
165\r
166 const Bool_t kDoESDFilter = kTRUE; //need this for JETAN with input ESDs\r
167 //const Bool_t kDoESDFilter = kFALSE;\r
168 if(kInputData == "ESD" && kDoESDFilter){\r
169 printf("Applying ESD filter cuts appropriate for jet analysis\n");\r
170 //\r
171 // Set of cuts\r
172 // \r
173 // standard\r
174 AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Loose");\r
175 //esdTrackCutsL->SetMinNClustersTPC(50);\r
176 //esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5);\r
b21a4af6 177 //esdTrackCutsL->SetRequireTPCRefit(kTRUE);\r
ef2c838f 178 //esdTrackCutsL->SetMaxDCAToVertexXY(2.4);\r
179 //esdTrackCutsL->SetMaxDCAToVertexZ(3.2);\r
180 //esdTrackCutsL->SetDCAToVertex2D(kTRUE);\r
181 //esdTrackCutsL->SetRequireSigmaToVertex(kFALSE);\r
b21a4af6 182 //esdTrackCutsL->SetAcceptKinkDaughters(kFALSE);\r
183 //\r
184 // hard\r
185 AliESDtrackCuts* esdTrackCutsH = new AliESDtrackCuts("AliESDtrackCuts", "Hard");\r
186 esdTrackCutsH->SetMinNClustersTPC(100);\r
187 esdTrackCutsH->SetMaxChi2PerClusterTPC(2.0);\r
b21a4af6 188 esdTrackCutsH->SetRequireTPCRefit(kTRUE);\r
ef2c838f 189 esdTrackCutsH->SetMaxDCAToVertexXY(2.4);\r
190 esdTrackCutsH->SetMaxDCAToVertexZ(3.2);\r
191 esdTrackCutsH->SetDCAToVertex2D(kTRUE);\r
192 esdTrackCutsH->SetRequireSigmaToVertex(kFALSE);\r
b21a4af6 193 esdTrackCutsH->SetAcceptKinkDaughters(kFALSE);\r
194 //\r
195 AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");\r
196 trackFilter->AddCuts(esdTrackCutsL);\r
197 // trackFilter->AddCuts(esdTrackCutsH);\r
198 //\r
199 AliAnalysisTaskESDfilter *esdfilter = new AliAnalysisTaskESDfilter("ESD Filter");\r
200 esdfilter->SetTrackFilter(trackFilter);\r
201 esdfilter->SetDebugLevel(10);\r
202 mgr->AddTask(esdfilter);\r
203 }\r
204\r
205\r
206 //-------------------------------------------------------------------------\r
207 //Define task, put here any other task that you want to use.\r
208 //-------------------------------------------------------------------------\r
209\r
210 //\r
211 // Jet analysis\r
212 //\r
d4df5eeb 213 if( kDoJetTask ){\r
214 AliAnalysisTaskJets *jetana = new AliAnalysisTaskJets("JetAnalysis",chain);\r
215 jetana->SetDebugLevel(2);\r
216 jetana->SetConfigFile(sconfig2); //Default ConfigJetAnalysisFastJet\r
217 //Uncommenting the following line produces too many AddAtAndExpand warnings for now.\r
218 if(kMergeAODs)jetana->ReadAODFromOutput(); //Uncomment when AOD merging\r
219 mgr->AddTask(jetana);\r
220 }\r
b21a4af6 221\r
222 //\r
223 // electron analysis\r
224 //\r
225 AliAnalysisTaskParticleCorrelation * taskpwg4 = new AliAnalysisTaskParticleCorrelation ("Particle");\r
226 taskpwg4->SetConfigFileName("ConfigAnalysisElectron"); //Default name is ConfigAnalysisElectron\r
227 mgr->AddTask(taskpwg4);\r
228 \r
229 // Create containers for input/output\r
230 AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();\r
231 AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();\r
232 AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("histos", TList::Class(),\r
233 AliAnalysisManager::kOutputContainer, "histos.root");\r
234\r
235\r
236 if(kInputData == "ESD" && kDoESDFilter){\r
237 mgr->ConnectInput (esdfilter, 0, cinput1 );\r
238 mgr->ConnectOutput (esdfilter, 0, coutput1 );\r
239 }\r
240\r
d4df5eeb 241 if( kDoJetTask ){\r
242 mgr->ConnectInput (jetana, 0, cinput1 );\r
243 mgr->ConnectOutput (jetana, 0, coutput1 );\r
244 mgr->ConnectOutput (jetana, 1, coutput2 );\r
245 }\r
246\r
b21a4af6 247 mgr->ConnectInput (taskpwg4, 0, cinput1 );\r
248 mgr->ConnectOutput (taskpwg4, 0, coutput1 );\r
249 mgr->ConnectOutput (taskpwg4, 1, coutput2 );\r
250 \r
251\r
252 //------------------------ \r
253 //Scaling task\r
254 //-----------------------\r
ef2c838f 255 cout<<">>> Scaling Task"<<endl;\r
256 Int_t nfiles = chain->GetListOfFiles()->GetEntriesFast();//chainxs->GetEntries();\r
b21a4af6 257 Int_t nevents = chain->GetEntries();\r
258 cout<<"Get? "<<kGetXSectionFromFileAndScale<<" nfiles "<<nfiles<<" nevents "<<nevents<<endl;\r
259 if(kGetXSectionFromFileAndScale && nfiles > 0){\r
ef2c838f 260 cout<<"Init AnaScale"<<endl;\r
b21a4af6 261 AliAnaScale * scale = new AliAnaScale("scale") ;\r
ef2c838f 262 cout<<"Summed xsection "<<xsection<<" Summed ntrials "<<ntrials<<" total events "<<nevents<<endl;\r
263 //Calculate the average\r
264 xsection/=nfiles;\r
265 ntrials/=nfiles;\r
266 cout<<"Average xsection "<<xsection<<" Average ntrials "<<ntrials<<" total events "<<nevents<<endl;\r
267 Double_t scaleFactor = xsection/ntrials/nevents ;\r
268 cout<<"Scale factor "<<scaleFactor<<endl;\r
269 scale->Set(scaleFactor) ;\r
b21a4af6 270 scale->MakeSumw2(kTRUE);//If you want histograms with error bars set to kTRUE\r
271 scale->SetDebugLevel(2);\r
272 mgr->AddTask(scale);\r
273 \r
ef2c838f 274 AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("histosscaled", \r
275 TList::Class(), AliAnalysisManager::kOutputContainer, "histosscaled.root");\r
b21a4af6 276 mgr->ConnectInput (scale, 0, coutput2);\r
277 mgr->ConnectOutput (scale, 0, coutput3 );\r
278 }\r
ef2c838f 279 \r
b21a4af6 280 //-----------------------\r
281 // Run the analysis\r
282 //----------------------- \r
283 TString smode = "";\r
284 if (mode==mLocal || mode == mLocalCAF) \r
285 smode = "local";\r
286 else if (mode==mPROOF) \r
287 smode = "proof";\r
288 else if (mode==mGRID) \r
289 smode = "local";\r
290 else if (mode==mPLUGIN) \r
291 smode = "grid";\r
292 \r
293 //mgr->ResetAnalysis();\r
294 mgr->InitAnalysis();\r
295 mgr->PrintStatus();\r
296 if (mode==mPROOF)\r
297 mgr->StartAnalysis(smode.Data(),kmydataset,1500,0);\r
298 else if (mode==mPLUGIN)\r
299 mgr->StartAnalysis(smode.Data());\r
300 else\r
ef2c838f 301 mgr->StartAnalysis(smode.Data(),chain);\r
b21a4af6 302\r
303 cout <<" Analysis ended sucessfully "<< endl ;\r
304\r
305 }\r
306 else cout << "Chain was not produced ! "<<endl;\r
307 \r
308}\r
309\r
310void LoadLibraries(const anaModes mode) {\r
311 \r
312 \r
313 //----------------------------------------------------------\r
314 // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< \r
315 //----------------------------------------------------------\r
316 if (mode==mLocal || mode == mLocalCAF || mode == mGRID || mode == mPLUGIN) {\r
b21a4af6 317\r
318 //--------------------------------------\r
319 // Load the needed libraries most of them already loaded by aliroot\r
320 //--------------------------------------\r
321 gSystem->Load("libTree.so");\r
322 gSystem->Load("libGeom.so");\r
323 gSystem->Load("libVMC.so");\r
324 gSystem->Load("libXMLIO.so");\r
d4df5eeb 325 if( kDoJetTask ){\r
326 gSystem->Load("libCGAL");\r
327 gSystem->Load("libfastjet");\r
328 gSystem->Load("libSISConePlugin");\r
329 }\r
b21a4af6 330\r
d4df5eeb 331 if(kUsePAR){\r
b21a4af6 332 //--------------------------------------------------------\r
333 //If you want to use root and par files from aliroot\r
334 //-------------------------------------------------------- \r
335 SetupPar("STEERBase");\r
336 SetupPar("ESD");\r
337 SetupPar("AOD");\r
338 SetupPar("ANALYSIS");\r
339 SetupPar("ANALYSISalice");\r
d4df5eeb 340 if( kDoJetTask ){\r
341 cerr<<"Now Loading JETAN"<<endl;\r
342 SetupPar("JETAN");\r
343 cerr<<"Done Loading JETAN"<<endl;\r
344 cerr<<"Now Loading FASTJETAN"<<endl;\r
345 SetupPar("FASTJETAN");\r
346 cerr<<"Done Loading FASTJETAN"<<endl;\r
347 }\r
b21a4af6 348 SetupPar("PWG4PartCorrBase");\r
349 SetupPar("PWG4PartCorrDep");\r
350 }\r
351 else{\r
352 //--------------------------------------------------------\r
353 // If you want to use already compiled libraries \r
354 // in the aliroot distribution\r
355 //--------------------------------------------------------\r
356 gSystem->Load("libSTEERBase");\r
357 gSystem->Load("libESD");\r
358 gSystem->Load("libAOD");\r
359 gSystem->Load("libANALYSIS");\r
360 gSystem->Load("libANALYSISalice");\r
d4df5eeb 361 if( kDoJetTask ){\r
362 gSystem->Load("libJETAN");\r
363 gSystem->Load("libFASTJETAN");\r
364 }\r
b21a4af6 365 gSystem->Load("libPWG4PartCorrBase");\r
366 gSystem->Load("libPWG4PartCorrDep");\r
367 }\r
368\r
369 \r
370 }\r
371\r
372 //---------------------------------------------------------\r
373 // <<<<<<<<<< PROOF mode >>>>>>>>>>>>\r
374 //---------------------------------------------------------\r
375 else if (mode==mPROOF) {\r
376 //\r
377 // Connect to proof\r
378 // Put appropriate username here\r
379 // char* myproofname = "alicecaf";\r
380 char* myproofname = "kread@localhost";\r
381\r
382 //TProof::Reset("proof://kread@lxb6046.cern.ch");\r
383 //TProof::Reset("proof://myproofname);\r
384 //TProof::Reset("myproofname",kTRUE);\r
385 gEnv->SetValue("XSec.GSI.DelegProxy","2"); \r
386 //TProof::Mgr(myproofname)->ShowROOTVersions();\r
387 //TProof::Mgr(myproofname)->SetROOTVersion("v5-23-04");\r
388 TProof::Open(myproofname);\r
389\r
390 // gProof->ClearPackages();\r
391 // gProof->SetLogLevel(5);\r
392 // gProof->ClearPackage("STEERBase");\r
393 // gProof->ClearPackage("ESD");\r
394 // gProof->ClearPackage("AOD");\r
395 // gProof->ClearPackage("ANALYSIS");\r
396 // gProof->ClearPackage("ANALYSISalice");\r
d4df5eeb 397 // if( kDoJetTask ){\r
398 // gProof->ClearPackage("JETAN");\r
399 // gProof->ClearPackage("FASTJETAN");\r
400 // }\r
b21a4af6 401 // gProof->ShowEnabledPackages();\r
402\r
403 // Enable the STEERBase Package\r
404 gProof->UploadPackage("STEERBase.par");\r
405 gProof->EnablePackage("STEERBase");\r
406 // Enable the ESD Package\r
407 gProof->UploadPackage("ESD.par");\r
408 gProof->EnablePackage("ESD");\r
409 // Enable the AOD Package\r
410 gProof->UploadPackage("AOD.par");\r
411 gProof->EnablePackage("AOD");\r
412 // Enable the Analysis Package\r
413 gProof->UploadPackage("ANALYSIS.par");\r
414 gProof->EnablePackage("ANALYSIS");\r
415 // Enable the Analysis Package\r
416 gProof->UploadPackage("ANALYSISalice.par");\r
417 gProof->EnablePackage("ANALYSISalice");\r
d4df5eeb 418 if( kDoJetTask ){\r
419 // Enable JETAN analysis\r
420 gProof->UploadPackage("JETAN.par");\r
421 gProof->EnablePackage("JETAN");\r
422 // Enable FASTJETAN analysis\r
423 gProof->UploadPackage("FASTJETAN.par");\r
424 gProof->EnablePackage("FASTJETAN");\r
425 }\r
b21a4af6 426\r
427 gProof->ShowEnabledPackages();\r
428 } \r
429 \r
430}\r
431\r
432void SetupPar(char* pararchivename)\r
433{\r
434 //Load par files, create analysis libraries\r
435 //For testing, if par file already decompressed and modified\r
436 //classes then do not decompress.\r
437 \r
438 TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; \r
439 TString parpar(Form("%s.par", pararchivename)) ; \r
440 if ( gSystem->AccessPathName(pararchivename) ) { \r
441 TString processline = Form(".! tar xvzf %s",parpar.Data()) ;\r
442 gROOT->ProcessLine(processline.Data());\r
443 }\r
444 \r
445 TString ocwd = gSystem->WorkingDirectory();\r
446 gSystem->ChangeDirectory(pararchivename);\r
447 \r
448 // check for BUILD.sh and execute\r
449 if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {\r
450 printf("*******************************\n");\r
451 printf("*** Building PAR archive ***\n");\r
452 cout<<pararchivename<<endl;\r
453 printf("*******************************\n");\r
454 \r
455 if (gSystem->Exec("PROOF-INF/BUILD.sh")) {\r
456 Error("runProcess","Cannot Build the PAR Archive! - Abort!");\r
457 return -1;\r
458 }\r
459 }\r
460 // check for SETUP.C and execute\r
461 if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {\r
462 printf("*******************************\n");\r
463 printf("*** Setup PAR archive ***\n");\r
464 cout<<pararchivename<<endl;\r
465 printf("*******************************\n");\r
466 gROOT->Macro("PROOF-INF/SETUP.C");\r
467 }\r
468 \r
469 gSystem->ChangeDirectory(ocwd.Data());\r
470 printf("Current dir: %s\n", ocwd.Data());\r
471}\r
472\r
473\r
474\r
475void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){\r
476 //Fills chain with data\r
ef2c838f 477\r
478 TString datafile="";\r
479 if(kInputData == "ESD") datafile = "AliESDs.root" ;\r
480 else if(kInputData == "AOD") datafile = "AliAOD.root" ;\r
481 else if(kInputData == "MC") datafile = "galice.root" ;\r
482\r
483 if(kInputData == "AOD") kXSFileName = "pyxsec_hists.root";\r
484\r
485 char * kGener = gSystem->Getenv("GENER");\r
486 if(kGener) {\r
487 cout<<"GENER "<<kGener<<endl;\r
488 if(!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";\r
489 else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";\r
490 else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;\r
491 }\r
492 else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;\r
493\r
494\r
b21a4af6 495 TString ocwd = gSystem->WorkingDirectory();\r
496 \r
497 //-----------------------------------------------------------\r
498 //Analysis of CAF data locally\r
499 //-----------------------------------------------------------\r
500 if(mode == mLocalCAF){\r
501 // Read the input list of files and add them to the chain\r
502 TString line;\r
503 ifstream in;\r
504 in.open("ESDlist.txt");\r
505 while (in.good()) {\r
506 in >> line;\r
507 if (line.Length() == 0) continue;\r
508 // cout << " line = " << line << endl;\r
509 chain->Add(line);\r
510 }\r
511 }\r
512 \r
513 //---------------------------------------\r
514 //Local files analysis\r
515 //---------------------------------------\r
516 else if(mode == mLocal){\r
517 //If you want to add several ESD files sitting in a common directory INDIR\r
518 //Specify as environmental variables the directory (INDIR), the number of files \r
519 //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)\r
520\r
521 if(gSystem->Getenv("INDIR")) \r
522 kInDir = gSystem->Getenv("INDIR") ; \r
523 else cout<<"INDIR not set, use default: "<<kInDir<<endl; \r
524 \r
525 if(gSystem->Getenv("PATTERN")) \r
526 kPattern = gSystem->Getenv("PATTERN") ; \r
527 else cout<<"PATTERN not set, use default: "<<kPattern<<endl;\r
528 \r
529 if(gSystem->Getenv("NFILES"))\r
530 kFile = atoi(gSystem->Getenv("NFILES")) ;\r
531 else cout<<"NFILES not set, use default: "<<kFile<<endl;\r
532 \r
533 //Check if env variables are set and are correct\r
534 if ( kInDir && kFile) {\r
535 printf("Get %d files from directory %s\n",kFile,kInDir);\r
536 if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist\r
537 printf("%s does not exist\n", kInDir) ;\r
538 return ;\r
539 }\r
540\r
541 //if(gSystem->Getenv("XSFILE")) \r
542 //kXSFileName = gSystem->Getenv("XSFILE") ; \r
543 //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl; \r
b21a4af6 544\r
545 cout<<"INDIR : "<<kInDir<<endl;\r
546 cout<<"NFILES : "<<kFile<<endl;\r
547 cout<<"PATTERN: " <<kPattern<<endl;\r
548 cout<<"XSFILE : "<<kXSFileName<<endl;\r
549 \r
ef2c838f 550\r
b21a4af6 551 //Loop on ESD files, add them to chain\r
552 Int_t event =0;\r
553 Int_t skipped=0 ; \r
554 char file[120] ;\r
555 char filexs[120] ;\r
556 \r
557 for (event = 0 ; event < kFile ; event++) {\r
558 sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; \r
559 sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ; \r
ef2c838f 560 TFile * dataFile = 0 ; \r
b21a4af6 561 //Check if file exists and add it, if not skip it\r
ef2c838f 562 if ( dataFile = TFile::Open(file)) {\r
563 if ( dataFile->Get(kTreeName) ) { \r
564 Int_t nEventsPerFile = ((TTree*) dataFile->Get(kTreeName)) ->GetEntries();\r
565 printf("++++ Adding %s, with %d events \n", file, nEventsPerFile) ;\r
b21a4af6 566 chain->AddFile(file);\r
ef2c838f 567 if(kGetXSectionFromFileAndScale) GetXsection(nEventsPerFile, filexs); \r
b21a4af6 568 }\r
569 }\r
570 else { \r
571 printf("---- Skipping %s\n", file) ;\r
572 skipped++ ;\r
573 }\r
574 }\r
575 printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ; \r
576 }\r
577 else {\r
578 TString input = "AliESDs.root" ;\r
579 cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;\r
580 chain->AddFile(input);\r
581 }\r
582 \r
583 }// local files analysis\r
584 \r
585 //------------------------------\r
586 //GRID xml files\r
587 //-----------------------------\r
588 else if(mode == mGRID){\r
589 //Get colection file. It is specified by the environmental\r
590 //variable XML\r
591\r
592 if(gSystem->Getenv("XML") )\r
593 kXML = gSystem->Getenv("XML");\r
594 else\r
595 sprintf(kXML, "collection.xml") ; \r
596 \r
597 if (!TFile::Open(kXML)) {\r
598 printf("No collection file with name -- %s -- was found\n",kXML);\r
599 return ;\r
600 }\r
601 else cout<<"XML file "<<kXML<<endl;\r
602\r
603 //Load necessary libraries and connect to the GRID\r
604 gSystem->Load("libNetx.so") ; \r
605 gSystem->Load("libRAliEn.so"); \r
606 TGrid::Connect("alien://") ;\r
607\r
608 //Feed Grid with collection file\r
609 //TGridCollection * collection = (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));\r
610 TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);\r
611 if (! collection) {\r
612 AliError(Form("%s not found", kXML)) ; \r
613 return kFALSE ; \r
614 }\r
615 TGridResult* result = collection->GetGridResult("",0 ,0);\r
616 \r
617 // Makes the ESD chain \r
618 printf("*** Getting the Chain ***\n");\r
ef2c838f 619 Int_t nEventsPerFile = 0;\r
b21a4af6 620 for (Int_t index = 0; index < result->GetEntries(); index++) {\r
621 TString alienURL = result->GetKey(index, "turl") ; \r
622 cout << "================== " << alienURL << endl ; \r
623 chain->Add(alienURL) ; \r
ef2c838f 624\r
625 if(kGetXSectionFromFileAndScale){\r
626 //Get the number of events per file.\r
627 //Do it only once, no need to open all the files.\r
628 if(i ==0 ) {\r
629 TFile * df = TFile::Open(alienURL);\r
630 nEventsPerFile = ((TTree*) df->Get(kTreeName)) ->GetEntries();\r
631 dataFile->Close();\r
632 } \r
633 alienURL.ReplaceAll(datafile,kXSFileName);\r
634 GetXsection(nEventsPerFile, alienURL);//chainxs->Add(alienURL) ; \r
635 }\r
b21a4af6 636 }\r
637 }// xml analysis\r
638\r
639 gSystem->ChangeDirectory(ocwd.Data());\r
640}\r
641\r
642//________________________________________________\r
ef2c838f 643void GetXsection(Int_t nEventsPerFile, TString filexs)\r
b21a4af6 644{\r
d4df5eeb 645 // Get the cross section from the corresponding file in the directory\r
ef2c838f 646 // where the data sits.\r
647 // The xsection and ntrials global variables are updated per each file.\r
648 // The average of these cuantities should be calculated after.\r
649\r
650 TFile *fxs = new TFile(filexs,"R");\r
651 if(kInputData =="AOD") { //needs improvement, in case of train with ESDs, reading output AODs for example this is wrong.\r
652 TList *l = (TList*) fxs->Get("cFilterList");\r
653 TH1F * hxs = l->FindObject("h1Xsec") ;\r
654 TH1F * htrial = l->FindObject("h1Trials") ;\r
655 if(htrial->GetEntries()!=hxs->GetEntries() || htrial->GetEntries()==0){\r
656 cout<<"Careful!!! Entries in histo for cross section "<<hxs->GetEntries()<< ", for trials "<<htrial->GetEntries()<<endl;\r
657 continue;\r
658 } \r
659 xsection += hxs->GetBinContent(1);\r
660 ntrials += htrial->GetBinContent(1)/htrial->GetEntries()/nEventsPerFile;\r
661 cout << "Chain: xsection " <<hxs->GetBinContent(1)<<" ntrials "<<htrial->GetBinContent(1)/htrial->GetEntries()<<endl; \r
662 }\r
663 else {\r
664 Double_t xs = 0;\r
665 UInt_t ntr = 0;\r
666 TTree * xstree = (TTree*)fxs->Get("Xsection");\r
667 xstree->SetBranchAddress("xsection",&xs);\r
668 xstree->SetBranchAddress("ntrials",&ntr);\r
669 xstree->GetEntry(0);\r
670 cout << "Chain: xsection " <<xs<<" ntrials "<<ntr<<endl;\r
671 xsection += xs ;\r
672 ntrials += ntr/nEventsPerFile; \r
673 } \r
b21a4af6 674}\r
675\r
676void ProcessEnvironment(){\r
677\r
d4df5eeb 678 if (gSystem->Getenv("anaInputData"))\r
679 kInputData = gSystem->Getenv("anaInputData");\r
680\r
ef2c838f 681 if (gSystem->Getenv("MODE"))\r
682 mode = atoi(gSystem->Getenv("MODE"));\r
b21a4af6 683\r
ef2c838f 684 if (gSystem->Getenv("CONFIG1"))\r
685 sprintf(sconfig1,gSystem->Getenv("CONFIG1"));\r
b21a4af6 686\r
ef2c838f 687 if (gSystem->Getenv("CONFIG2"))\r
688 sprintf(sconfig2,gSystem->Getenv("CONFIG2"));\r
b21a4af6 689\r
ef2c838f 690 if (gSystem->Getenv("CONFIG3"))\r
691 sprintf(sconfig3,gSystem->Getenv("CONFIG3"));\r
b21a4af6 692\r
ef2c838f 693 if (gSystem->Getenv("SEVENT"))\r
694 sevent = atoi (gSystem->Getenv("SEVENT"));\r
b21a4af6 695 \r
d4df5eeb 696 printf("Process: Variables: kInputData %s, mode %d, config2 %s, config3 %s, sevent %d\n",kInputData.Data(),mode,sconfig2,sconfig3,sevent);\r
b21a4af6 697\r
b21a4af6 698}\r