has to be positive otherwise no clusters made
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / macros / emcalReclusterize.C
1 /* $Id:  $ */
2 //--------------------------------------------------
3 // Example macro to do Calorimeters filtering
4 // copy ESDs into AODs
5 //
6 // Pay attention to the options and definitions
7 // set in the lines below
8 //
9 //  Author : Gustavo Conesa Balbastre (INFN-LNF)
10 //
11 //-------------------------------------------------
12 enum anaModes {mLocal, mLocalCAF,mPROOF,mGRID};
13 //mLocal: Analyze locally files in your computer
14 //mLocalCAF: Analyze locally CAF files
15 //mPROOF: Analyze CAF files with PROOF
16
17 //---------------------------------------------------------------------------
18 //Settings to read locally several files, only for "mLocal" mode
19 //The different values are default, they can be set with environmental 
20 //variables: INDIR, PATTERN, NFILES, respectivelly
21 char * kInDir = "/user/data/files/"; 
22 char * kPattern = ""; // Data are in files kInDir/kPattern+i 
23 Int_t kFile = 1; // Number of files
24 //---------------------------------------------------------------------------
25 //Collection file for grid analysis
26 char * kXML = "collection.xml";
27 //---------------------------------------------------------------------------
28
29 const TString kInputData = "ESD"; //ESD, AOD, MC
30 TString kTreeName = "esdTree";
31 Bool_t kUsePhysSel = kFALSE;
32
33 void emcalReclusterize(Int_t mode=mLocal)
34 {
35   // Main
36   char cmd[200] ; 
37   sprintf(cmd, ".! rm -rf aod.root") ; 
38   gROOT->ProcessLine(cmd) ; 
39   //--------------------------------------------------------------------
40   // Load analysis libraries
41   // Look at the method below, 
42   // change whatever you need for your analysis case
43   // ------------------------------------------------------------------
44   LoadLibraries(mode) ;
45   //gSystem->Unload("libPWG4CaloCalib.so");
46   //Try to set the new library
47   //gSystem->Load("./PWG4CaloCalib/libPWG4CaloCalib.so");
48   //gSystem->ListLibraries();
49
50   //-------------------------------------------------------------------------------------------------
51   //Create chain from ESD and from cross sections files, look below for options.
52   //-------------------------------------------------------------------------------------------------
53   if(kInputData == "ESD") kTreeName = "esdTree" ;
54   else if(kInputData == "AOD") kTreeName = "aodTree" ;
55   else {
56     cout<<"Wrong  data type "<<kInputData<<endl;
57     break;
58   }
59
60   TChain *chain       = new TChain(kTreeName) ;
61   CreateChain(mode, chain);  
62
63   if(chain){
64     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
65     
66     //--------------------------------------
67     // Make the analysis manager
68     //-------------------------------------
69     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
70
71     // AOD output handler
72     AliAODHandler* aodoutHandler   = new AliAODHandler();
73     aodoutHandler->SetOutputFileName("aod.root");
74     ////aodoutHandler->SetCreateNonStandardAOD();
75     mgr->SetOutputEventHandler(aodoutHandler);
76     
77     //input
78     if(kInputData == "ESD")
79       {
80         // ESD handler
81         AliESDInputHandler *esdHandler = new AliESDInputHandler();
82         mgr->SetInputEventHandler(esdHandler);
83         esdHandler->SetReadFriends(kFALSE);
84         cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
85       }
86     if(kInputData == "AOD")
87       {
88         // AOD handler
89         AliAODInputHandler *aodHandler = new AliAODInputHandler();
90         mgr->SetInputEventHandler(aodHandler);
91         cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
92         
93       }
94     
95     mgr->SetDebugLevel(1);
96     
97     //-------------------------------------------------------------------------
98     //Define task, put here any other task that you want to use.
99     //-------------------------------------------------------------------------
100     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
101     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
102     
103     // ESD physics selection task
104     if(kInputData == "ESD" && kUsePhysSel){
105       gROOT->LoadMacro("AddTaskPhysicsSelection.C");
106       AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
107     }
108     
109     // Create containers for input/output
110     AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
111     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
112  
113     AliAnalysisTaskEMCALClusterize * clusterize = new AliAnalysisTaskEMCALClusterize("EMCALClusterize");
114     if(kUsePhysSel)clusterize->SelectCollisionCandidates();
115     //clusterize->SetOCDBPath("local://$ALICE_ROOT/OCDB");
116     clusterize->FillAODFile(kFALSE); // fill aod.root with clusters?, not really needed for analysis.
117     clusterize->JustUnfold(kTRUE); // if TRUE, do just unfolding, do not recluster cells
118     AliEMCALRecParam * params = clusterize->GetRecParam();
119     params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN);
120     params->SetClusteringThreshold(0.1); // 100 MeV                                             
121     params->SetMinECut(0.01);  //10 MeV    
122     params->SetUnfold(kFALSE);
123     params->SetW0(4.5);
124     params->SetTimeCut(1e6);///Open this cut for AODs
125     params->SetTimeMin(-1);//Open this cut for AODs
126     params->SetTimeMax(1e6);//Open this cut for AODs
127
128 //    TGeoHMatrix *matrix[4];
129     
130 //     double rotationMatrix[4][9] = {-0.014585, -0.999892, -0.002031, 0.999892, -0.014589,  0.001950, -0.001979, -0.002003,  0.999996,
131 //                                 -0.014585,  0.999892,  0.002031, 0.999892,  0.014589, -0.001950, -0.001979,  0.002003, -0.999996,
132 //                                 -0.345861, -0.938280, -0.003412, 0.938281, -0.345869,  0.001950, -0.003010, -0.002527,  0.999992,
133 //                                 -0.345861,  0.938280,  0.003412, 0.938281,  0.345869, -0.001950, -0.003010,  0.002527, -0.999992};
134     
135 //     double translationMatrix[4][3] = {0.367264,    446.508738,  175.97185+0.3,
136 //                                    1.078181,    445.826258, -174.026758+0.3,
137 //                                    -153.843916, 418.304256,  175.956905+0.8,
138 //                                    -152.649580, 417.621779, -174.040392+0.8};
139
140
141 //     double rotationMatrix[4][9] = {-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996,
142 //                               -0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996,
143 //                               -0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990,
144 //                               -0.345864,  0.938278,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990};
145     
146 //     double translationMatrix[4][3] = {0.351659,    447.576446,  176.269742,
147 //                                    1.062577,    446.893974, -173.728870,
148 //                                    -154.213287, 419.306156,  176.753692,
149 //                                    -153.018950, 418.623681, -173.243605};
150
151 //     for(int j=0; j<4; j++)
152 //       {
153 //      matrix[j] = new TGeoHMatrix();
154 //      matrix[j]->SetRotation(rotationMatrix[j]);
155 //      matrix[j]->SetTranslation(translationMatrix[j]);
156 //      matrix[j]->Print();
157 //      clusterize->SetGeometryMatrixInSM(matrix[j],j);
158 //       }
159     
160     
161 //     clusterize->SwitchOnLoadOwnEMCALGeometryMatrices();
162     
163         
164     mgr->AddTask(clusterize);
165     mgr->ConnectInput  (clusterize,  0, cinput1);
166     mgr->ConnectOutput (clusterize, 0, coutput1 );
167    
168
169     //-----------------------
170     // Run the analysis
171     //-----------------------    
172     TString smode = "";
173     if (mode==mLocal || mode == mLocalCAF) 
174       smode = "local";
175     else if (mode==mPROOF) 
176       smode = "proof";
177     else if (mode==mGRID) 
178       smode = "local";
179     
180     mgr->InitAnalysis();
181     mgr->PrintStatus();
182     mgr->StartAnalysis(smode.Data(),chain);
183
184 cout <<" Analysis ended sucessfully "<< endl ;
185
186   }
187   else cout << "Chain was not produced ! "<<endl;
188   
189 }
190
191 void  LoadLibraries(const anaModes mode) {
192   
193   //--------------------------------------
194   // Load the needed libraries most of them already loaded by aliroot
195   //--------------------------------------
196   gSystem->Load("libTree.so");
197   gSystem->Load("libGeom.so");
198   gSystem->Load("libVMC.so");
199   gSystem->Load("libXMLIO.so");
200   gSystem->Load("libMatrix.so");
201   gSystem->Load("libPhysics.so");  
202
203   //----------------------------------------------------------
204   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
205   //----------------------------------------------------------
206   if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
207     //--------------------------------------------------------
208     // If you want to use already compiled libraries 
209     // in the aliroot distribution
210     //--------------------------------------------------------
211
212     gSystem->Load("libSTEERBase.so");
213     gSystem->Load("libESD.so");
214     gSystem->Load("libAOD.so");
215     gSystem->Load("libANALYSIS.so");
216     gSystem->Load("libANALYSISalice.so");
217     gSystem->Load("libANALYSISalice.so");
218     gSystem->Load("libEMCALUtils.so");  
219
220     //SetupPar("EMCALUtils"); 
221     //SetupPar("PWG4CaloCalib"); 
222     
223     TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
224     gSystem->Load("libPWG4CaloCalib.so");
225    
226     /*
227       //     gSystem->Load("libPWG4omega3pi.so");
228       //     gSystem->Load("libCORRFW.so");
229       //     gSystem->Load("libPWG3base.so");
230       //     gSystem->Load("libPWG3muon.so");
231  */
232     //--------------------------------------------------------
233     //If you want to use root and par files from aliroot
234     //--------------------------------------------------------  
235     /*     
236            SetupPar("STEERBase");
237            SetupPar("ESD");
238            SetupPar("AOD");
239            SetupPar("ANALYSIS");
240            SetupPar("ANALYSISalice");  
241            SetupPar("PHOSUtils");
242            SetupPar("EMCALUtils");
243            //Create Geometry
244            TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
245            SetupPar("PWG4CaloCalib");
246 */
247   }
248
249   //---------------------------------------------------------
250   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
251   //---------------------------------------------------------
252   else if (mode==mPROOF) {
253     //
254     // Connect to proof
255     // Put appropriate username here
256     // TProof::Reset("proof://mgheata@lxb6046.cern.ch"); 
257     TProof::Open("proof://mgheata@lxb6046.cern.ch");
258     
259     //    gProof->ClearPackages();
260     //    gProof->ClearPackage("ESD");
261     //    gProof->ClearPackage("AOD");
262     //    gProof->ClearPackage("ANALYSIS");   
263     //    gProof->ClearPackage("PWG4PartCorrBase");
264     //    gProof->ClearPackage("PWG4PartCorrDep");
265     
266     // Enable the STEERBase Package
267     gProof->UploadPackage("STEERBase.par");
268     gProof->EnablePackage("STEERBase");
269     // Enable the ESD Package
270     gProof->UploadPackage("ESD.par");
271     gProof->EnablePackage("ESD");
272     // Enable the AOD Package
273     gProof->UploadPackage("AOD.par");
274     gProof->EnablePackage("AOD");
275     // Enable the Analysis Package
276     gProof->UploadPackage("ANALYSIS.par");
277     gProof->EnablePackage("ANALYSIS");
278     // Enable the PHOS geometry Package
279     //gProof->UploadPackage("PHOSUtils.par");
280     //gProof->EnablePackage("PHOSUtils");
281     // Enable PartCorr analysis
282     gProof->UploadPackage("PWG4PartCorrBase.par");
283     gProof->EnablePackage("PWG4PartCorrBase");
284     gProof->UploadPackage("PWG4PartCorrDep.par");
285     gProof->EnablePackage("PWG4PartCorrDep");    
286     gProof->ShowEnabledPackages();
287   }  
288   
289 }
290
291 void SetupPar(char* pararchivename)
292 {
293   //Load par files, create analysis libraries
294   //For testing, if par file already decompressed and modified
295   //classes then do not decompress.
296  
297   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
298   TString parpar(Form("%s.par", pararchivename)) ; 
299   if ( gSystem->AccessPathName(parpar.Data()) ) {
300     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
301     TString processline(Form(".! make %s", parpar.Data())) ; 
302     gROOT->ProcessLine(processline.Data()) ;
303     gSystem->ChangeDirectory(cdir) ; 
304     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
305     gROOT->ProcessLine(processline.Data()) ;
306   } 
307   if ( gSystem->AccessPathName(pararchivename) ) {  
308     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
309     gROOT->ProcessLine(processline.Data());
310   }
311   
312   TString ocwd = gSystem->WorkingDirectory();
313   gSystem->ChangeDirectory(pararchivename);
314   
315   // check for BUILD.sh and execute
316   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
317     printf("*******************************\n");
318     printf("*** Building PAR archive    ***\n");
319     cout<<pararchivename<<endl;
320     printf("*******************************\n");
321     
322     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
323       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
324       return -1;
325     }
326   }
327   // check for SETUP.C and execute
328   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
329     printf("*******************************\n");
330     printf("*** Setup PAR archive       ***\n");
331     cout<<pararchivename<<endl;
332     printf("*******************************\n");
333     gROOT->Macro("PROOF-INF/SETUP.C");
334   }
335   
336   gSystem->ChangeDirectory(ocwd.Data());
337   printf("Current dir: %s\n", ocwd.Data());
338 }
339
340
341
342 void CreateChain(const anaModes mode, TChain * chain){
343   //Fills chain with data
344   TString ocwd = gSystem->WorkingDirectory();
345   
346   //-----------------------------------------------------------
347   //Analysis of CAF data locally and with PROOF
348   //-----------------------------------------------------------
349   if(mode ==mPROOF || mode ==mLocalCAF){
350     // Chain from CAF
351     gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
352     // The second parameter is the number of input files in the chain
353     chain = CreateESDChain("ESD12001.txt", 5);  
354   }
355   
356   //---------------------------------------
357   //Local files analysis
358   //---------------------------------------
359   else if(mode == mLocal){    
360     //If you want to add several ESD files sitting in a common directory INDIR
361     //Specify as environmental variables the directory (INDIR), the number of files 
362     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
363
364     if(gSystem->Getenv("INDIR"))  
365       kInDir = gSystem->Getenv("INDIR") ; 
366     else cout<<"INDIR not set, use default: "<<kInDir<<endl;    
367     
368     if(gSystem->Getenv("PATTERN"))   
369       kPattern = gSystem->Getenv("PATTERN") ; 
370     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
371     
372     if(gSystem->Getenv("NFILES"))
373       kFile = atoi(gSystem->Getenv("NFILES")) ;
374     else cout<<"NFILES not set, use default: "<<kFile<<endl;
375     
376     //Check if env variables are set and are correct
377     if ( kInDir  && kFile) {
378       printf("Get %d files from directory %s\n",kFile,kInDir);
379       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
380         printf("%s does not exist\n", kInDir) ;
381         return ;
382       }
383
384   
385       cout<<"INDIR   : "<<kInDir<<endl;
386       cout<<"NFILES  : "<<kFile<<endl;
387       cout<<"PATTERN : " <<kPattern<<endl;
388
389       TString datafile="";
390       if(kInputData == "ESD") datafile = "AliESDs.root" ;
391       else if(kInputData == "AOD") datafile = "AliAODs.root" ;
392       else if(kInputData == "MC")  datafile = "galice.root" ;
393       
394       //Loop on ESD files, add them to chain
395       Int_t event =0;
396       Int_t skipped=0 ; 
397       char file[120] ;
398       
399       for (event = 0 ; event < kFile ; event++) {
400         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; 
401         TFile * fESD = 0 ; 
402         //Check if file exists and add it, if not skip it
403         if ( fESD = TFile::Open(file)) {
404           if ( fESD->Get(kTreeName) ) { 
405             printf("++++ Adding %s\n", file) ;
406             chain->AddFile(file);
407           }
408         }
409         else { 
410           printf("---- Skipping %s\n", file) ;
411           skipped++ ;
412         }
413       }
414       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;      
415     }
416     else {
417       TString input = "AliESDs.root" ;
418       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
419       chain->AddFile(input);
420     }
421     
422   }// local files analysis
423   
424   //------------------------------
425   //GRID xml files
426   //-----------------------------
427   else if(mode == mGRID){
428     //Get colection file. It is specified by the environmental
429     //variable XML
430
431     if(gSystem->Getenv("XML") )
432       kXML = gSystem->Getenv("XML");
433     else
434       sprintf(kXML, "collection.xml") ; 
435     
436     if (!TFile::Open(kXML)) {
437       printf("No collection file with name -- %s -- was found\n",kXML);
438       return ;
439     }
440     else cout<<"XML file "<<kXML<<endl;
441
442     //Load necessary libraries and connect to the GRID
443     gSystem->Load("libNetx.so") ; 
444     gSystem->Load("libRAliEn.so"); 
445     TGrid::Connect("alien://") ;
446
447     //Feed Grid with collection file
448     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
449     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
450     if (! collection) {
451       AliError(Form("%s not found", kXML)) ; 
452       return kFALSE ; 
453     }
454     TGridResult* result = collection->GetGridResult("",0 ,0);
455    
456     // Makes the ESD chain 
457     printf("*** Getting the Chain       ***\n");
458     for (Int_t index = 0; index < result->GetEntries(); index++) {
459       TString alienURL = result->GetKey(index, "turl") ; 
460       cout << "================== " << alienURL << endl ; 
461       chain->Add(alienURL) ; 
462     }
463   }// xml analysis
464   
465   gSystem->ChangeDirectory(ocwd.Data());
466 }
467