db54c49c713b2ac494aabf3517b6067b6abe2363
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / macros / ana.C
1 /* $Id:  $ */
2 //--------------------------------------------------
3 // Example macro to do analysis with the 
4 // analysis classes in CaloTrackCorrelations
5 // Can be executed with Root and AliRoot
6 //
7 // Pay attention to the options and definitions
8 // set in the lines below
9 //
10 //  Author : Gustavo Conesa Balbastre (INFN-LNF)
11 //
12 //-------------------------------------------------
13 enum anaModes {mLocal=0, mPROOF=1, mPlugin=2, mGRID=3};
14 //mLocal    = 0: Analyze locally files in your computer
15 //mPROOF    = 1: Analyze files on GRID with Plugin
16 //mPlugin   = 2: Analyze files on GRID with Plugin
17 //mGRID     = 3: Analyze files on GRID, jobs launched from aliensh
18
19 //---------------------------------------------------------------------------
20 // Settings to read locally several files, only for "mLocal" mode
21 // The different values are default, they can be set with environmental 
22 // variables: INDIR, PATTERN, NFILES, respectivelly
23
24 char * kInDir   = "/user/data/files/"; 
25 char * kPattern = ""; // Data are in files kInDir/kPattern+i 
26 Int_t  kFile    = 1; 
27
28 //---------------------------------------------------------------------------
29 // Dataset for proof analysis, mode=mPROOF
30 // char * kDataset = "/alice/vernet/PbPb_LHC10h_ESD";
31
32 char *  kDatasetPROOF     = "/alice/vernet/LHC11b_149646";
33 //char *  kDatasetPROOF     = "/alice/vernet/LHC11b10a_AOD046";//LHC11d_AOD076
34 Int_t   kDatasetNMaxFiles = 20;
35 TString ccin2p3UserName   = "arbor" ;
36 TString alienUserName     = "narbor" ;
37
38 //---------------------------------------------------------------------------
39 // Collection file for grid analysis
40
41 char * kXML = "collection.xml";
42
43 //---------------------------------------------------------------------------
44 //Scale histograms from file. Change to kTRUE when xsection file exists
45 //Put name of file containing xsection 
46 //Put number of events per ESD file
47 //This is an specific case for normalization of Pythia files.
48 const Bool_t kGetXSectionFromFileAndScale = kFALSE ;
49 const char * kXSFileName = "pyxsec.root";
50
51 //---------------------------------------------------------------------------
52
53 //Set some default values, but used values are set in the code!
54
55 Bool_t  kMC        = kFALSE; //With real data kMC = kFALSE
56 TString kInputData = "ESD"; //ESD, AOD, MC, deltaAOD
57 Int_t   kYear      = 2011;
58 TString kCollision = "pp";
59 Bool_t  outAOD     = kTRUE; //Some tasks doesnt need it.
60 TString kTreeName;
61 TString kPass      = "";
62 char    kTrigger[1024];
63 Int_t   kRun       = 0;
64
65 //________________________
66 void ana(Int_t mode=mGRID)
67 {
68   // Main
69   
70   //--------------------------------------------------------------------
71   // Load analysis libraries
72   // Look at the method below, 
73   // change whatever you need for your analysis case
74   // ------------------------------------------------------------------
75   
76   LoadLibraries(mode) ;
77   //gSystem->ListLibraries();
78   
79   //-------------------------------------------------------------------------------------------------
80   //Create chain from ESD and from cross sections files, look below for options.
81   //-------------------------------------------------------------------------------------------------
82   
83   // Set kInputData and kTreeName looking to the kINDIR
84   
85   CheckInputData(mode);
86   
87   // Check global analysis settings  
88   
89   CheckEnvironmentVariables();
90   
91   printf("*********************************************\n");
92   printf("*** Input data < %s >, pass %s, tree < %s >, MC?  < %d > ***\n",kInputData.Data(),kPass.Data(),kTreeName.Data(),kMC);
93   printf("*********************************************\n");
94   
95   TChain * chain   = new TChain(kTreeName) ;
96   TChain * chainxs = new TChain("Xsection") ;
97   CreateChain(mode, chain, chainxs); 
98   
99   printf("*********************************************\n");
100   printf("number of entries # %lld, skipped %d\n", chain->GetEntries()) ;       
101   printf("*********************************************\n");
102   
103   if(!chain)
104   { 
105     printf("STOP, no chain available\n"); 
106     return;
107   }
108   
109   AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
110   
111   //------------------------------------------
112   //  Alien handler part
113   //------------------------------------------
114   AliAnalysisGrid *alienHandler=0x0;
115   if(mode==mPlugin)
116   {
117     // Create and configure the alien handler plugin
118     gROOT->LoadMacro("CreateAlienHandler.C");
119     alienHandler = CreateAlienHandler();
120     if (!alienHandler) return;
121   }  
122   
123   //--------------------------------------
124   // Make the analysis manager
125   //-------------------------------------
126   AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
127   //AliAnalysisManager::SetUseProgressBar(kTRUE);
128   //mgr->SetSkipTerminate(kTRUE);
129   //mgr->SetNSysInfo(1);
130   
131   if(mode==mPlugin)
132   {
133     // Connect plugin to the analysis manager
134     mgr->SetGridHandler(alienHandler);
135   }
136   
137   // MC handler
138   if((kMC || kInputData == "MC") && !kInputData.Contains("AOD"))
139   {
140     AliMCEventHandler* mcHandler = new AliMCEventHandler();
141     mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
142     mgr->SetMCtruthEventHandler(mcHandler);
143     if( kInputData == "MC") 
144     {
145       cout<<"MC INPUT EVENT HANDLER"<<endl;
146       mgr->SetInputEventHandler(NULL);
147     }
148   }
149   
150   
151   // AOD output handler
152   if(kInputData!="deltaAOD" && outAOD)
153   {
154     cout<<"Init output handler"<<endl;
155     AliAODHandler* aodoutHandler   = new AliAODHandler();
156     aodoutHandler->SetOutputFileName("aod.root");
157     ////aodoutHandler->SetCreateNonStandardAOD();
158     mgr->SetOutputEventHandler(aodoutHandler);
159   }
160   
161   //input
162   
163   if(kInputData == "ESD")
164   {
165     // ESD handler
166     AliESDInputHandler *esdHandler = new AliESDInputHandler();
167     esdHandler->SetReadFriends(kFALSE);
168     mgr->SetInputEventHandler(esdHandler);
169     cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
170   }
171   else if(kInputData.Contains("AOD"))
172   {
173     // AOD handler
174     AliAODInputHandler *aodHandler = new AliAODInputHandler();
175     mgr->SetInputEventHandler(aodHandler);
176     if(kInputData == "deltaAOD") aodHandler->AddFriend("deltaAODCaloTrackCorr.root");
177     cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
178   }
179   //mgr->RegisterExternalFile("deltaAODCaloTrackCorr.root");
180   //mgr->SetDebugLevel(1); // For debugging, do not uncomment if you want no messages.
181   
182   TString outputFile = AliAnalysisManager::GetCommonFileName(); 
183   
184   //-------------------------------------------------------------------------
185   //Define task, put here any other task that you want to use.
186   //-------------------------------------------------------------------------
187   
188   // Physics selection
189   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); 
190   AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kMC); 
191   
192   // Centrality
193   if(kCollision=="PbPb")
194   {
195     gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
196     AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
197   }
198   
199   // Simple event counting tasks
200   AddTaskCounter("");   // All
201   //AddTaskCounter("MB"); // Min Bias
202   if(!kMC)
203   {
204     AddTaskCounter("INT7"); // Min Bias
205     //AddTaskCounter("EMC1"); // Trig Th > 1.5 GeV approx
206     AddTaskCounter("EMC7"); // Trig Th > 4-5 GeV 
207     AddTaskCounter("PHOS"); //  
208   }
209     
210   // -----------------
211   // Photon conversion
212   // ----------------- 
213 /*  
214   if(kInputData=="ESD"){
215     printf("* Configure photon conversion analysis in macro \n");
216     TString arguments = "-run-on-train -use-own-xyz  -force-aod -mc-off ";
217     gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/ConfigGammaConversion.C");
218     AliAnalysisTaskGammaConversion * taskGammaConversion = 
219     ConfigGammaConversion(arguments,mgr->GetCommonInputContainer());
220     taskGammaConversion->SelectCollisionCandidates();
221     
222     // Gamma Conversion AOD to AODPWG4Particle
223     AliAnalysisTaskGCPartToPWG4Part * taskGCToPC = new AliAnalysisTaskGCPartToPWG4Part("GCPartToPWG4Part");
224     taskGCToPC->SetGammaCutId("90035620401003321022000000090");
225     mgr->AddTask(taskGCToPC);
226     mgr->ConnectInput  (taskGCToPC, 0, mgr->GetCommonInputContainer() );
227     mgr->ConnectOutput (taskGCToPC, 0, mgr->GetCommonOutputContainer()); 
228   }
229 */  
230   
231   Bool_t kPrint   = kFALSE;
232   Bool_t deltaAOD = kFALSE;
233   gROOT->LoadMacro("AddTaskCaloTrackCorr.C");   // $ALICE_ROOT/PWG4/macros
234   gROOT->LoadMacro("AddTaskEMCALClusterize.C"); // $ALICE_ROOT/PWG4/CaloCalib/macros  
235   
236   
237   // ------
238   // Tracks
239   // ------  
240
241   // Track isolation-correlation analysis and EMCAL QA analysis
242   AliAnalysisTaskCaloTrackCorrelation *anamb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
243                                                                kYear,kRun,kCollision,"INT7","");   // PHOS trigger
244   
245   AliAnalysisTaskCaloTrackCorrelation *anatr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
246                                                                kYear,kRun,kCollision,"EMC7","");  
247   
248  
249   // -----
250   // EMCAL
251   // -----  
252   
253   Bool_t  bTrackMatch = kTRUE;
254   Int_t   minEcell    = 50;  // 50  MeV (10 MeV used in reconstruction)
255   Int_t   minEseed    = 100; // 100 MeV
256   Int_t   dTime       = 0;   // default, 250 ns
257   Int_t   wTime       = 0;   // default 425 < T < 825 ns
258   TString clTrigger   = "";  // Do not select, do Min Bias and triggered
259   
260   //Analysis with clusterizer V1
261   AliAnalysisTaskEMCALClusterize * clv1 = AddTaskEMCALClusterize(kMC,"V1",clTrigger,kRun,kPass, bTrackMatch,
262                                                                  minEcell,minEseed,dTime,wTime);    
263   
264   TString arrayNameV1(Form("V1_Ecell%d_Eseed%d_DT%d_WT%d",minEcell,minEseed, dTime,wTime));
265   printf("Name of clusterizer array: %s\n",arrayNameV1.Data());
266   
267   if(!kMC)
268   {
269     
270     AliAnalysisTaskCaloTrackCorrelation *anav1tr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
271                                                                    kYear,kRun,kCollision,"EMC7",arrayNameV1);
272     
273     AliAnalysisTaskCaloTrackCorrelation *anav1mb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
274                                                                    kYear,kRun,kCollision,"INT7",arrayNameV1);
275   }
276   else 
277   {// No trigger (should be MB, but for single particle productions it does not work)
278     
279     AliAnalysisTaskCaloTrackCorrelation *anav1  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
280                                                                  kYear,kRun,kCollision,"",arrayNameV1);
281   }
282   
283   
284   
285   //Analysis with clusterizer V2
286   AliAnalysisTaskEMCALClusterize * clv2 = AddTaskEMCALClusterize(kMC,"V2",clTrigger,kRun,kPass, bTrackMatch,
287                                                                  minEcell,minEseed,dTime,wTime);    
288   
289   TString arrayNameV2(Form("V2_Ecell%d_Eseed%d_DT%d_WT%d",minEcell,minEseed, dTime,wTime));
290   printf("Name of clusterizer array: %s\n",arrayNameV2.Data());
291   
292   if(!kMC)
293   {
294     
295     AliAnalysisTaskCaloTrackCorrelation *anav2tr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
296                                                                    kYear,kRun,kCollision,"EMC7",arrayNameV2);
297     
298     AliAnalysisTaskCaloTrackCorrelation *anav2mb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
299                                                                    kYear,kRun,kCollision,"INT7",arrayNameV2);
300   }
301   else 
302   {// No trigger (should be MB, but for single particle productions it does not work)
303     AliAnalysisTaskCaloTrackCorrelation *anav2  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
304                                                                  kYear,kRun,kCollision,"",arrayNameV2);    
305   }
306  
307   /*
308   // -----
309   // PHOS
310   // -----
311   
312   //Add here PHOS tender or whatever is needed
313   
314   if(!kMC)
315   {
316     
317     AliAnalysisTaskCaloTrackCorrelation *anav1tr = AddTaskCaloTrackCorr(kInputData, "PHOS", kPrint,kMC, deltaAOD,  outputFile.Data(), 
318                                                                   kYear,kRun,kCollision,"PHOS",""); // PHOS trigger
319     
320     AliAnalysisTaskCaloTrackCorrelation *anav1mb = AddTaskCaloTrackCorr(kInputData, "PHOS",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
321                                                                   kYear,kRun,kCollision,"INT7","");
322     
323   }
324   else 
325   {// No trigger
326     
327     AliAnalysisTaskCaloTrackCorrelation *anav1mb = AddTaskCaloTrackCorr(kInputData, "PHOS",   kPrint,kMC, deltaAOD,  outputFile.Data(), 
328                                                                   kYear,kRun,kCollision,"","");
329     
330   }
331
332   */
333   //-----------------------
334   // Run the analysis
335   //-----------------------    
336   mgr->InitAnalysis();
337   mgr->PrintStatus();
338   
339   if      (mode == mPlugin) mgr->StartAnalysis("grid");
340   else if (mode == mPROOF ) mgr->StartAnalysis("proof",chain);
341   else                      mgr->StartAnalysis("local",chain);
342   
343   cout <<" Analysis ended sucessfully "<< endl ;
344   
345 }
346
347 //_____________________________
348 void  LoadLibraries(Int_t mode)
349 {
350   
351   if (mode == mPROOF) {
352     //TProof::Mgr("ccalpmaster")->SetROOTVersion("ALICE_v5-27-06b");
353     gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/EnableAliRootForLAF.C");
354     TProof* proof = EnableAliRootForLAF("ccaplmaster",nPROOFWorkers.Data(),ccin2p3UserName.Data(),alienUserName.Data(),"",kFALSE,kTRUE,kTRUE,"OADB:ANALYSIS:ANALYSISalice:AOD:ESD:CORRFW:STEERBase:EMCALUtils:PHOSUtils:PWG4PartCorrBase:PWG4PartCorrDep:PWG4CaloCalib");
355     
356     //  TProof* proof = TProof::Open("ccaplmaster",Form("workers=%s",nPROOFWorkers.Data()));
357     
358     //     //proof->ClearPackages();
359     //     proof->UploadPackage("STEERBase");
360     //     proof->UploadPackage("ESD");
361     //     proof->UploadPackage("AOD");
362     //     proof->UploadPackage("ANALYSIS");
363     //     proof->UploadPackage("OADB");
364     //     proof->UploadPackage("ANALYSISalice");
365     //     proof->UploadPackage("CORRFW");
366     //     //proof->UploadPackage("JETAN");
367     //     proof->UploadPackage("PHOSUtils");
368     //     proof->UploadPackage("EMCALUtils");
369     //     proof->UploadPackage("PWG4PartCorrBase");
370     //     proof->UploadPackage("PWG4PartCorrDep");
371     //     proof->UploadPackage("PWG4CaloCalib");
372     
373     //     proof->EnablePackage("STEERBase");
374     //     proof->EnablePackage("ESD");
375     //     proof->EnablePackage("AOD");
376     //     proof->EnablePackage("ANALYSIS");
377     //     proof->EnablePackage("OADB");
378     //     proof->EnablePackage("ANALYSISalice");
379     //     proof->EnablePackage("CORRFW");
380     //     //proof->EnablePackage("JETAN");
381     //     proof->EnablePackage("PHOSUtils");
382     //     proof->EnablePackage("EMCALUtils");
383     //     proof->EnablePackage("PWG4PartCorrBase");
384     //     proof->EnablePackage("PWG4PartCorrDep");
385     //     proof->EnablePackage("PWG4CaloCalib");
386     return;
387   }  
388   
389   //--------------------------------------
390   // Load the needed libraries most of them already loaded by aliroot
391   //--------------------------------------
392   gSystem->Load("libTree.so");
393   gSystem->Load("libGeom.so");
394   gSystem->Load("libVMC.so");
395   gSystem->Load("libXMLIO.so");
396   gSystem->Load("libMatrix.so");
397   gSystem->Load("libPhysics.so");
398   gSystem->Load("libMinuit.so"); // Root + libraries to if reclusterization is done
399   
400   gSystem->Load("libSTEERBase.so");
401   gSystem->Load("libGui.so"); // Root + libraries to if reclusterization is done
402   gSystem->Load("libCDB.so"); // Root + libraries to if reclusterization is done
403   gSystem->Load("libESD.so"); // Root + libraries to if reclusterization is done
404   gSystem->Load("libAOD.so");
405   gSystem->Load("libRAWDatabase.so"); // Root + libraries to if reclusterization is done
406   gSystem->Load("libProof.so"); 
407   gSystem->Load("libANALYSIS.so");
408   gSystem->Load("libSTEER.so"); // Root + libraries to if reclusterization is done
409   
410   gSystem->Load("libRAWDatarec.so"); // Root + libraries to if reclusterization is done
411   gSystem->Load("libRAWDatasim.so"); // Root + libraries to if reclusterization is done
412   gSystem->Load("libVZERObase.so");  // Root + libraries to if reclusterization is done
413   gSystem->Load("libVZEROrec.so");   // Root + libraries to if reclusterization is done
414   
415   gSystem->Load("libEMCALUtils");
416   //SetupPar("EMCALUtils");
417   gSystem->Load("libEMCALraw");  // Root + libraries to if reclusterization is done
418   gSystem->Load("libEMCALbase"); // Root + libraries to if reclusterization is done
419   gSystem->Load("libEMCALsim");  // Root + libraries to if reclusterization is done
420   gSystem->Load("libEMCALrec");  // Root + libraries to if reclusterization is done
421   //SetupPar("EMCALraw");
422   //SetupPar("EMCALbase");
423   //SetupPar("EMCALsim");
424   //SetupPar("EMCALrec");
425   
426   gSystem->Load("libANALYSISalice.so");
427   //gSystem->Load("libTENDER.so"); 
428   //gSystem->Load("libTENDERSupplies.so");
429   
430   gSystem->Load("libPHOSUtils");
431   gSystem->Load("libEMCALUtils");
432   gSystem->Load("libPWG4PartCorrBase");
433   gSystem->Load("libPWG4PartCorrDep");
434   gSystem->Load("libPWG4CaloCalib");
435   //SetupPar("PWG4PartCorrBase");
436   //SetupPar("PWG4PartCorrDep");
437   //SetupPar("PWG4CaloCalib");
438   
439   //gSystem->Load("libJETAN");
440   //gSystem->Load("FASTJETAN");
441   //gSystem->Load("PWG4JetTasks");
442
443   gSystem->Load("libCORRFW.so");
444   gSystem->Load("libPWG4GammaConv.so"); 
445   //SetupPar("PWG4GammaConv"); 
446   
447   // needed for plugin?
448   gSystem->AddIncludePath("-I$ALICE_ROOT");
449   gSystem->AddIncludePath("-I./");     
450   
451 }
452
453 //_________________________________
454 void SetupPar(char* pararchivename)
455 {
456   //Load par files, create analysis libraries
457   //For testing, if par file already decompressed and modified
458   //classes then do not decompress.
459   
460   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
461   TString parpar(Form("%s.par", pararchivename)) ; 
462   
463   if ( gSystem->AccessPathName(pararchivename) ) {  
464     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
465     gROOT->ProcessLine(processline.Data());
466   }
467   
468   TString ocwd = gSystem->WorkingDirectory();
469   gSystem->ChangeDirectory(pararchivename);
470   
471   // check for BUILD.sh and execute
472   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
473     printf("*******************************\n");
474     printf("*** Building PAR archive    ***\n");
475     cout<<pararchivename<<endl;
476     printf("*******************************\n");
477     
478     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
479       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
480       return -1;
481     }
482   }
483   // check for SETUP.C and execute
484   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
485     printf("*******************************\n");
486     printf("*** Setup PAR archive       ***\n");
487     cout<<pararchivename<<endl;
488     printf("*******************************\n");
489     gROOT->Macro("PROOF-INF/SETUP.C");
490   }
491   
492   gSystem->ChangeDirectory(ocwd.Data());
493   printf("Current dir: %s\n", ocwd.Data());
494 }
495
496 //______________________________________
497 void CheckInputData(const anaModes mode)
498 {
499   //Sets input data and tree
500   
501   TString ocwd = gSystem->WorkingDirectory();
502   
503   //---------------------------------------
504   //Local files analysis
505   //---------------------------------------
506   if(mode == mLocal){    
507     //If you want to add several ESD files sitting in a common directory INDIR
508     //Specify as environmental variables the directory (INDIR), the number of files 
509     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
510     
511     if(gSystem->Getenv("INDIR"))  
512       kInDir = gSystem->Getenv("INDIR") ; 
513     else cout<<"INDIR not set, use default: "<<kInDir<<endl;    
514     
515     TString sindir(kInDir);
516     if     (sindir.Contains("pass1")) kPass = "pass1";
517     else if(sindir.Contains("pass2")) kPass = "pass2";
518     else if(sindir.Contains("pass3")) kPass = "pass3";
519     
520     if(gSystem->Getenv("PATTERN"))   
521       kPattern = gSystem->Getenv("PATTERN") ; 
522     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
523     
524     cout<<"INDIR   : "<<kInDir<<endl;
525     cout<<"NFILES  : "<<kFile<<endl;
526     
527     char fileE[120] ;   
528     char fileA[120] ;   
529     char fileG[120] ;
530     char fileEm[120] ;   
531     for (Int_t event = 0 ; event < kFile ; event++) {
532       sprintf(fileE,  "%s/%s%d/AliESDs.root",    kInDir,kPattern,event) ; 
533       sprintf(fileA,  "%s/%s%d/AliAOD.root",     kInDir,kPattern,event) ; 
534       sprintf(fileG,  "%s/%s%d/galice.root",     kInDir,kPattern,event) ; 
535       sprintf(fileEm, "%s/%s%d/embededAOD.root", kInDir,kPattern,event) ; 
536       
537       TFile * fESD = TFile::Open(fileE) ; 
538       TFile * fAOD = TFile::Open(fileA) ; 
539       
540       //Check if file exists and add it, if not skip it
541       if (fESD) 
542       {
543         kTreeName  = "esdTree";
544         kInputData = "ESD";
545         TFile * fG = TFile::Open(fileG);
546         if(fG) { kMC = kTRUE; fG->Close();}
547         else     kMC = kFALSE;
548         
549         // Get run number
550         TTree* esdTree = (TTree*)fESD->Get("esdTree");
551         AliESDEvent* esd = new AliESDEvent();
552         esd->ReadFromTree(esdTree);
553         esdTree->GetEvent(0);
554         kRun = esd->GetRunNumber();
555         
556         return;
557       }
558       else if(fAOD)
559       {
560         kTreeName  = "aodTree";
561         kInputData = "AOD";
562         if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
563         else kMC = kFALSE;
564         
565         // Get run number
566         TTree* aodTree = (TTree*)fAOD->Get("aodTree");
567         AliAODEvent* aod = new AliAODEvent();
568         aod->ReadFromTree(aodTree);
569         aodTree->GetEvent(0);
570         kRun = aod->GetRunNumber();
571         return;
572       }
573       else if(TFile::Open(fileEm))
574       {
575         kTreeName  = "aodTree";
576         kInputData = "AOD";
577         kMC        = kTRUE;
578         
579         return;
580       }
581       else if(TFile::Open(fileG))
582       {
583         kTreeName  = "TE";
584         kInputData = "MC";
585         kMC        = kTRUE;
586         return;
587       }
588     }
589     
590     if(fESD) fESD->Close();
591     if(fAOD) fAOD->Close();
592     
593   }// local files analysis
594   
595   //------------------------------
596   //GRID xml files
597   //-----------------------------
598   else if(mode == mGRID){
599     //Get colection file. It is specified by the environmental
600     //variable XML
601     
602     if(gSystem->Getenv("XML") )
603       kXML = gSystem->Getenv("XML");
604     else
605       sprintf(kXML, "collection.xml") ; 
606     
607     if (!TFile::Open(kXML)) {
608       printf("No collection file with name -- %s -- was found\n",kXML);
609       return ;
610     }
611     else cout<<"XML file "<<kXML<<endl;
612     
613     //Load necessary libraries and connect to the GRID
614     gSystem->Load("libNetx.so") ; 
615     gSystem->Load("libRAliEn.so"); 
616     TGrid::Connect("alien://") ;
617     
618     //Feed Grid with collection file
619     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
620     if (! collection) {
621       AliError(Form("%s not found", kXML)) ; 
622       return kFALSE ; 
623     }
624     TGridResult* result = collection->GetGridResult("",0 ,0);
625     
626     for (Int_t index = 0; index < result->GetEntries(); index++) {
627       TString alienURL = result->GetKey(index, "turl") ; 
628       cout << "================== " << alienURL << endl ; 
629       
630       if     (alienURL.Contains("pass1")) kPass = "pass1";
631       else if(alienURL.Contains("pass2")) kPass = "pass2";
632       else if(alienURL.Contains("pass3")) kPass = "pass3";
633       
634       kRun = AliAnalysisManager::GetRunFromAlienPath(alienURL.Data());
635       printf("Run number from alien path = %d\n",kRun);
636       
637       TFile * fAOD = 0 ; 
638       //Check if file exists and add it, if not skip it
639       if (alienURL.Contains("AliESDs.root"))  
640       {
641         kTreeName  = "esdTree";
642         kInputData = "ESD";
643         alienURL.ReplaceAll("AliESDs.root","galice.root");
644         if(TFile::Open(alienURL)) kMC=kTRUE;
645         else kMC = kFALSE;
646         return;
647       }
648       else if(alienURL.Contains("AliAOD.root"))
649       {
650         kTreeName  = "aodTree";
651         kInputData = "AOD";
652         fAOD = TFile::Open(alienURL);
653         if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
654         else kMC = kFALSE;
655         return;
656       }
657       else if(alienURL.Contains("embededAOD.root"))
658       {
659         kTreeName  = "aodTree";
660         kInputData = "AOD";
661         kMC=kTRUE;
662         return;
663       }
664       else if(alienURL.Contains("galice.root"))
665       {
666         kTreeName  = "TE";
667         kInputData = "MC";
668         kMC=kTRUE;
669         return;
670       } 
671     }
672   }// xml analysis
673   //------------------------------
674   //PROOF files
675   //-----------------------------
676   else if(mode == mPROOF){
677     
678     TFileCollection* coll  = gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
679     
680     TIter iter(coll->GetList());
681     
682     TFileInfo* fileInfo = 0;
683     while ((fileInfo = dynamic_cast<TFileInfo*> (iter())))
684     {
685       if (fileInfo->GetFirstUrl()) {
686         TString ProofURL = fileInfo->GetFirstUrl()->GetUrl();
687         cout << "================== " << ProofURL << endl ; 
688         
689         if     (ProofURL.Contains("pass1")) kPass = "pass1";
690         else if(ProofURL.Contains("pass2")) kPass = "pass2";
691         else if(ProofURL.Contains("pass3")) kPass = "pass3";
692         
693         kRun = AliAnalysisManager::GetRunFromAlienPath(ProofURL.Data());
694         printf("Run number from alien path = %d\n",kRun);
695         
696         TFile * fAOD = 0 ; 
697         //Check if file exists and add it, if not skip it
698         if (ProofURL.Contains("AliESDs.root"))  
699         {
700           kTreeName  = "esdTree";
701           kInputData = "ESD";
702           alienURL.ReplaceAll("AliESDs.root","galice.root");
703           if(TFile::Open(ProofURL)) kMC=kTRUE;
704           else kMC = kFALSE;
705           
706           return;
707         }
708         else if(ProofURL.Contains("AliAOD.root"))
709         {
710           kTreeName  = "aodTree";
711           kInputData = "AOD";
712           fAOD = TFile::Open(ProofURL);
713           if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
714           else kMC = kFALSE;
715           return;
716         }
717         else if(ProofURL.Contains("embededAOD.root"))
718         {
719           kTreeName  = "aodTree";
720           kInputData = "AOD";
721           kMC=kTRUE;
722           return;
723         }
724         else if(ProofURL.Contains("galice.root"))
725         {
726           kTreeName  = "TE";
727           kInputData = "MC";
728           kMC=kTRUE;
729           return;
730         } 
731       }
732     }
733   }// proof analysis
734   
735   gSystem->ChangeDirectory(ocwd.Data());
736   
737 }
738
739 //_____________________________________________________________________
740 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs)
741 {
742   //Fills chain with data
743   TString ocwd = gSystem->WorkingDirectory();
744   
745   //---------------------------------------
746   // Local files analysis
747   //---------------------------------------
748   if(mode == mLocal){    
749     //If you want to add several ESD files sitting in a common directory INDIR
750     //Specify as environmental variables the directory (INDIR), the number of files 
751     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
752     
753     if(gSystem->Getenv("INDIR"))  
754       kInDir = gSystem->Getenv("INDIR") ; 
755     else cout<<"INDIR not set, use default: "<<kInDir<<endl;    
756     
757     if(gSystem->Getenv("PATTERN"))   
758       kPattern = gSystem->Getenv("PATTERN") ; 
759     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
760     
761     if(gSystem->Getenv("NFILES"))
762       kFile = atoi(gSystem->Getenv("NFILES")) ;
763     else cout<<"NFILES not set, use default: "<<kFile<<endl;
764     
765     //Check if env variables are set and are correct
766     if ( kInDir  && kFile) {
767       printf("Get %d files from directory %s\n",kFile,kInDir);
768       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
769         printf("%s does not exist\n", kInDir) ;
770         return ;
771       }
772       
773       //if(gSystem->Getenv("XSFILE"))  
774       //kXSFileName = gSystem->Getenv("XSFILE") ; 
775       //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;   
776       char * kGener = gSystem->Getenv("GENER");
777       if(kGener) {
778         cout<<"GENER "<<kGener<<endl;
779         if     (!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";
780         else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";
781         else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;
782       }
783       else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;
784       
785       cout<<"INDIR   : "<<kInDir     <<endl;
786       cout<<"NFILES  : "<<kFile      <<endl;
787       cout<<"PATTERN : "<<kPattern   <<endl;
788       cout<<"XSFILE  : "<<kXSFileName<<endl;
789       
790       TString datafile="";
791       if     (kInputData == "ESD")        datafile = "AliESDs.root" ;
792       else if(kInputData.Contains("AOD")) datafile = "AliAOD.root"  ;
793       else if(kInputData == "MC")         datafile = "galice.root"  ;
794       
795       //Loop on ESD/AOD/MC files, add them to chain
796       Int_t event =0;
797       Int_t skipped=0 ; 
798       char file[120] ;
799       char filexs[120] ;
800       
801       for (event = 0 ; event < kFile ; event++) {
802         sprintf(file,   "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; 
803         sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ; 
804         TFile * fData = 0 ; 
805         //Check if file exists and add it, if not skip it
806         if ( fData = TFile::Open(file)) {
807           if ( fData->Get(kTreeName) ) { 
808             printf("++++ Adding %s\n", file) ;
809             chain->AddFile(file);
810             chainxs->Add(filexs) ; 
811           }
812         }
813         else { 
814           printf("---- Skipping %s\n", file) ;
815           skipped++ ;
816         }
817       }
818     }
819     else {
820       TString input = "AliESDs.root" ;
821       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
822       chain->AddFile(input);
823     }
824     
825   }// local files analysis
826   
827   //------------------------------
828   // GRID xml files
829   //------------------------------
830   else if(mode == mGRID){
831     //Get colection file. It is specified by the environmental
832     //variable XML
833     
834     //Feed Grid with collection file
835     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
836     if (! collection) {
837       AliError(Form("%s not found", kXML)) ; 
838       return kFALSE ; 
839     }
840     
841     TGridResult* result = collection->GetGridResult("",0 ,0);
842     
843     // Makes the ESD chain 
844     printf("*** Getting the Chain       ***\n");
845     for (Int_t index = 0; index < result->GetEntries(); index++) {
846       TString alienURL = result->GetKey(index, "turl") ; 
847       cout << "================== " << alienURL << endl ; 
848       chain->Add(alienURL) ; 
849       alienURL.ReplaceAll("AliESDs.root",kXSFileName);
850       chainxs->Add(alienURL) ; 
851     }
852   }// xml analysis
853   
854   //------------------------------
855   // PROOF
856   //------------------------------
857   else if (mode == mPROOF) {
858     
859     TFileCollection* ds= gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
860     
861     gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/dataset_management/CreateChainFromDataSet.C");
862     chain = CreateChainFromDataSet(ds, kTreeName , kDatasetNMaxFiles);
863     printf("chain has %d entries\n",chain->GetEntries());
864   }
865   
866   gSystem->ChangeDirectory(ocwd.Data());
867   
868 }
869
870 //_________________________________________________________________
871 void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
872 {
873   // Read the PYTHIA statistics from the file pyxsec.root created by
874   // the function WriteXsection():
875   // integrated cross section (xsection) and
876   // the  number of Pyevent() calls (ntrials)
877   // and calculate the weight per one event xsection/ntrials
878   // The spectrum calculated by a user should be
879   // multiplied by this weight, something like this:
880   // TH1F *userSpectrum ... // book and fill the spectrum
881   // userSpectrum->Scale(weight)
882   //
883   // Yuri Kharlov 19 June 2007
884   // Gustavo Conesa 15 April 2008
885   Double_t xsection = 0;
886   UInt_t    ntrials = 0;
887   xs = 0;
888   ntr = 0;
889   
890   Int_t nfiles =  tree->GetEntries()  ;
891   if (tree && nfiles > 0) 
892   {
893     tree->SetBranchAddress("xsection",&xsection);
894     tree->SetBranchAddress("ntrials",&ntrials);
895     for(Int_t i = 0; i < nfiles; i++){
896       tree->GetEntry(i);
897       xs += xsection ;
898       ntr += ntrials ;
899       cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl; 
900     }
901     
902     xs =   xs /  nfiles;
903     ntr =  ntr / nfiles;
904     cout << "-----------------------------------------------------------------"<<endl;
905     cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl; 
906     cout << "-----------------------------------------------------------------"<<endl;
907   } 
908   else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
909   
910 }
911
912 //______________________________
913 void CheckEnvironmentVariables()
914 {
915   
916   sprintf(kTrigger,"");
917   
918   Bool_t bRecalibrate = kFALSE;
919   Bool_t bBadChannel = kFALSE;
920   
921   for (int i=0; i< gApplication->Argc();i++){
922     
923 #ifdef VERBOSEARGS
924     
925     printf("Arg  %d:  %s\n",i,gApplication->Argv(i));
926     
927 #endif
928     
929     TString sRun = "";
930     
931     if (!(strcmp(gApplication->Argv(i),"--trigger")))
932       sprintf(trigger,gApplication->Argv(i+1));
933     
934     if (!(strcmp(gApplication->Argv(i),"--recalibrate")))
935       bRecalibrate = atoi(gApplication->Argv(i+1));
936     
937     if (!(strcmp(gApplication->Argv(i),"--badchannel")))
938       bBadChannel = atoi(gApplication->Argv(i+1));
939     
940     if (!(strcmp(gApplication->Argv(i),"--run"))){
941       sRun = gApplication->Argv(i+1);
942       if(sRun.Contains("LHC10")) {
943         kYear = 2010;
944       }
945       else {
946         if(kRun <=0){
947           kRun = atoi(gApplication->Argv(i+1));
948         }
949         else printf("** Run number already set  to %d, do not set to %d\n",kRun,atoi(gApplication->Argv(i+1)));
950       }//numeric run
951     }//--run available
952     
953   }// args loop
954   
955   if(!sRun.Contains("LHC10")){
956     if ( kRun < 140000) {
957       kYear = 2010;
958       if( kRun >= 136851 ) kCollision = "PbPb";
959     }
960     else{
961       kYear = 2011;
962     }
963   }
964   
965   if(kMC) sprintf(kTrigger,"");
966   
967   printf("*********************************************\n");
968   //printf("*** Settings trigger %s, recalibrate %d, remove bad channels %d, year %d, collision %s, run %d ***\n",
969   //       kTrigger,bRecalibrate,bBadChannel, kYear,kCollision.Data(), kRun);
970   printf("*** Settings year %d, collision %s, run %d ***\n",kYear,kCollision.Data(), kRun);
971   printf("*********************************************\n");
972   
973 }
974
975 //_______________________________________________
976 void AddTaskCounter(const TString trigger = "MB")
977 {
978   
979   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
980   
981   AliAnalysisTaskCounter * counter =  new AliAnalysisTaskCounter(Form("Counter%s",trigger.Data()));
982   if(kRun > 140000 && kRun < 146900) counter ->RejectFastCluster();
983   if     (kCollision=="pp"  )   counter->SetZVertexCut(50.);  //Open cut
984   else if(kCollision=="PbPb")   counter->SetZVertexCut(10.);  //Centrality defined in this range.
985   
986   if(trigger=="EMC7")
987   {
988     printf("counter trigger EMC7\n");
989     counter->SelectCollisionCandidates(AliVEvent::kEMC7);
990   }
991   else if (trigger=="INT7")
992   {
993     printf("counter trigger INT7\n");
994     counter->SelectCollisionCandidates(AliVEvent::kINT7);
995   }
996   if(trigger=="EMC1")
997   {
998     printf("counter trigger EMC1\n");
999     counter->SelectCollisionCandidates(AliVEvent::kEMC1);
1000   }
1001   else if(trigger=="MB")
1002   {
1003     printf("counter trigger MB\n");
1004     counter->SelectCollisionCandidates(AliVEvent::kMB);
1005   }
1006   else if(trigger=="PHOS")
1007   {
1008     printf("counter trigger PHOS\n");
1009     counter->SelectCollisionCandidates(AliVEvent::kPHI7);
1010   }
1011   
1012   TString outputFile = AliAnalysisManager::GetCommonFileName(); 
1013   AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
1014   
1015   AliAnalysisDataContainer *coutput = 
1016   mgr->CreateContainer(Form("Counter%s",trigger.Data()), TList::Class(), AliAnalysisManager::kOutputContainer,  outputFile.Data());
1017   mgr->AddTask(counter);
1018   mgr->ConnectInput  (counter, 0, cinput1);
1019   mgr->ConnectOutput (counter, 1, coutput);
1020   
1021 }
1022
1023
1024
1025
1026