New macros for analysis train ESDtoAOD, compliant with analysis framework, possibilit...
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 7 Jul 2009 14:08:31 +0000 (14:08 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 7 Jul 2009 14:08:31 +0000 (14:08 +0000)
PWG3/READMEmuon
PWG3/muon/AddTaskTagCreation.C [new file with mode: 0644]
PWG3/muon/AnalysisTrainFromESDToAOD.C [new file with mode: 0644]
PWG3/muon/CreateAlienHandler_FromESDToAOD.C [new file with mode: 0644]
PWG3/muon/ReadAOD_MCBranch.C [new file with mode: 0644]

index f126330..d781de5 100644 (file)
@@ -9,11 +9,38 @@ The following filter is used in the official analysis train, in order to copy th
 1) AliAnalysisTaskESDMuonFilter.h  --> analysis task to copy the muon information from the ESD to the standard AOD
 2) AliAnalysisTaskESDMuonFilter.cxx --> analysis task to copy the muon information from the ESD to the standard AOD
 
-In order to test locally the analysis train from the ESD to the Standard AOD production, the following macro is provided
-
-1) AnalysisTrainMuonLocal.C 
-   - The input file is the ESD
-   - The outputs are the standard AOD and the AOD tag file
+In order to test the analysis train from the ESD to the Standard AOD production, the following macro is provided:
+
+  AnalysisTrainFromESDToAOD.C
+
+Input file: ESD
+Output files: standard AOD (+AOD tags)
+
+Two wagons are attached to the analysis train: 
+   1) $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C 
+   2) $ALICE_ROOT/PWG3/muon/AddTaskTagCreation.C 
+
+Several flags can be activated: 
+   1) iESDfilter: to activate, using AddTaskESDFilter.C, the copy of the ESD information into the AOD.
+      In AddTaskESDFilter.C there is the possibility to apply cuts on the tracks and muon tracks in order
+      to reject them before filling the AOD. 
+   2) iAODAddMCBranch: to activate the inclusion of the MC branch (containing Kinematics info) into a branch of the AOD 
+   3) iAODTagCreation: to activate, using AddTaskTagCreation.C, the AOD tag creation
+      
+Runninng options tested: 
+   1) GRID (with/without AliEn plugin)
+   2) LOCAL 
+
+If the AliEn plugin is required, the macro 
+   CreateAlienHandler_FromESDToAOD.C
+can be used, in order to configure the jdl.
+
+AnalysisTrainFromESDToAOD.C is an updated version of the macro AnalysisTrainMuonLocal.C 
+(input file: ESD, output files standard AOD and the AOD tag file)
+The main differences are related to implementations to be compliant with the analysis framework,
+possibility to run on the grid with/without the Alien Plugin and possibility of adding the MC truth in a branch of the AOD.
+
+The macro ReadAOD_MCBranch.C is an example on how to access MC information stored in the AOD.
    
 ===================================================
   Creation of the MUON-AOD from the standard AOD
diff --git a/PWG3/muon/AddTaskTagCreation.C b/PWG3/muon/AddTaskTagCreation.C
new file mode 100644 (file)
index 0000000..beccac1
--- /dev/null
@@ -0,0 +1,38 @@
+AliAnalysisTaskTagCreator *AddTaskTagCreation()
+{
+
+// Creates tag AOD files
+
+   // Get the pointer to the existing analysis manager via the static access method.
+   //==============================================================================
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if (!mgr) {
+      ::Error("AddTaskTagCreator ", "No analysis manager to connect to.");
+      return NULL;
+   }   
+   
+   // Check input/output handlers
+   //===============================================================================
+   AliAODInputHandler *aod_input_h = (AliAODInputHandler*)mgr->GetInputEventHandler();
+   if (!aod_input_h) {
+      ::Error("AddTaskTagCreator ", "Input AOD handler does not exist!.");
+      return NULL;
+   }
+      
+   // Create the task, add it to the manager and configure it.
+   //===========================================================================   
+
+   AliAnalysisTaskTagCreator *tagcreator  = new AliAnalysisTaskTagCreator("TagCreator");
+   mgr->AddTask(tagcreator);
+   
+   // Create ONLY the output containers for the data produced by the task.
+   // Get and connect other common input/output containers via the manager as below
+   //==============================================================================
+    // Tag container                                                         
+    AliAnalysisDataContainer *cout_tags = mgr->CreateContainer("cTag",TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root");
+    
+    mgr->ConnectInput  (tagcreator ,  0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput (tagcreator ,  1, cout_tags);
+
+   return tagcreator ;
+}   
diff --git a/PWG3/muon/AnalysisTrainFromESDToAOD.C b/PWG3/muon/AnalysisTrainFromESDToAOD.C
new file mode 100644 (file)
index 0000000..0d3b1e8
--- /dev/null
@@ -0,0 +1,317 @@
+class AliAnalysisGrid; 
+const char *dataset   = "";
+TString anaLibs = "";
+Int_t iESDfilter       = 1;
+Int_t iAODTagCreation  = 1;
+Int_t iAODAddMCBranch  = 1;
+
+void AnalysisTrainFromESDToAOD(const char *analysisMode = "GRID", Bool_t usePLUGIN = kFALSE, Int_t nev=12345678)
+
+//========================================================================
+// The macro produces a standard AOD starting from ESD files. 
+// (Simplified version of ANALYSIS/macros/AnalysisTrainNew.C)
+//
+// Two wagons are attached to the train: 
+//      AddTaskESDFilter.C and AddTaskTagCreation.C 
+//
+// If the iESDfilter flag is activated, in AddTaskESDFilter.C 
+// two tasks are executed :
+// 1- with the first one (AliAnalysisTaskESDfilter), 
+//    all the branches of the AOD are filled apart from the muons 
+// 2- with the second task (AliAnalysisTaskESDMuonFilter) 
+//    muons tracks are added to the tracks branch 
+//
+// In AddTaskESDFilter.C there is the possibility to apply cuts 
+// on the tracks and muon tracks in order to reject them before 
+// filling the AOD. 
+//
+// - if the flag iAODAddMCBranch is activated the MC branch 
+// (containing Kinematics info) is added to the AOD 
+
+// - if the iAODTagCreation flag is activated, in AddTaskTagCreation.C the 
+//  AliAnalysisTaskTagCreator task is executed in order to create aod tags.  
+//
+//
+// Options tested: (case sensitive)
+//    GRID (with/without AliEn plugin)
+//    LOCAL (you have to specify in TChain *CreateChain(...) 
+//          the directory where your data are)
+//========================================================================
+
+{
+    // Global configuration flags 
+    //=====================================================================
+    Bool_t debug         = kTRUE;
+    Bool_t readTR        = kFALSE;      
+    Bool_t useKFILTER    = kFALSE;  // add MC Branch 
+    if(iAODAddMCBranch) useKFILTER=kTRUE;
+    if(strcmp(analysisMode,"LOCAL")==0) usePLUGIN = kFALSE; 
+       
+    // Load common libraries (STEERBase, ESD, AOD, ANALYSIS. ANALYSISalice)
+    //=====================================================================
+    LoadCommonLibraries(analysisMode);
+    gROOT->ProcessLine(".include $ALICE_ROOT/include");
+   
+    // Load analysis specific libraries
+    //=====================================================================
+    if (iESDfilter) {
+       if(!strcmp(analysisMode, "LOCAL")){
+         gSystem->Load("libPWG3base.so");
+         gSystem->Load("libPWG3muon.so");
+       } 
+       else if(!strcmp(analysisMode, "GRID"))LoadAnalysisLibraries(analysisMode);
+     }
+     
+    // If Plugin is used, load macro with JDL parameters
+    //=====================================================================
+    if(usePLUGIN){
+      gROOT->LoadMacro("CreateAlienHandler_FromESDToAOD.C");
+      AliAnalysisGrid *alienHandler = CreateAlienHandler_FromESDToAOD();  
+      if (!alienHandler) return;     
+    }
+
+     // Create the chain. This is dependent on the analysis mode.
+     //=====================================================================
+     if(!usePLUGIN){
+       if (!strcmp(analysisMode, "GRID")) TGrid::Connect("alien://");
+       TChain* chain = CreateChain(analysisMode,""); 
+     }
+     
+    // Create the train and set-up the handlers
+    //=====================================================================
+    AliAnalysisManager *mgr  = new AliAnalysisManager("Analysis Train", "Analysis Train for standard AOD production");
+
+    // GRID handler
+    if(usePLUGIN) mgr->SetGridHandler(alienHandler); 
+    
+    // ESD input handler
+       AliESDInputHandler *esdHandler = new AliESDInputHandler();
+       mgr->SetInputEventHandler(esdHandler);      
+
+    // Monte Carlo handler
+    if (iAODAddMCBranch) {
+       AliMCEventHandler* mcHandler = new AliMCEventHandler();
+       mgr->SetMCtruthEventHandler(mcHandler);
+       mcHandler->SetReadTR(readTR);
+    } 
+     
+    // AOD output handler
+    AliAODHandler* aodHandler   = new AliAODHandler();
+    mgr->SetOutputEventHandler(aodHandler);
+    aodHandler->SetOutputFileName("AliAODs.root");
+       
+    // Debugging if requested
+    if (debug) mgr->SetDebugLevel(3);
+       
+    // Load the tasks configuration macros for all included wagons
+    //=====================================================================    
+    if (iESDfilter) {
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
+       AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(useKFILTER);
+       
+       if (iAODTagCreation) {
+// use this line if AddTaskTagCreation.C is available in the grid aliroot version
+       if(!strcmp(analysisMode, "LOCAL")){
+         gROOT->LoadMacro("$ALICE_ROOT/PWG3/muon/AddTaskTagCreation.C");
+       } else  {
+// uncomment  this line if AddTaskTagCreation.C is available in the grid aliroot version
+//       gROOT->LoadMacro("$ALICE_ROOT/PWG3/muon/AddTaskTagCreation.C");
+// otherwise temporary: (and AddTaskTagCreation.C must be also added in the jdl)
+         gROOT->LoadMacro("AddTaskTagCreation.C");
+       }        
+       AliAnalysisTaskTagCreator *tagcreator = AddTaskTagCreation();
+      }       
+    }  
+
+    // Run the analysis
+    //=====================================================================    
+    if (mgr->InitAnalysis()) {
+       mgr->PrintStatus();
+       if(usePLUGIN) mgr->StartAnalysis("GRID");
+       else mgr->StartAnalysis("local", chain,nev);
+    }  
+}
+
+//______________________________________________________________________________
+Bool_t LoadCommonLibraries(const char *mode)
+{
+// Load common analysis libraries.
+   Int_t imode = -1;
+   if (!strcmp(mode, "LOCAL")) imode = 0;
+   if (!strcmp(mode, "PROOF")) imode = 1;
+   if (!strcmp(mode, "GRID"))  imode = 2;
+   if (!gSystem->Getenv("ALICE_ROOT")) {
+      ::Error("LoadCommonLibraries", "Analysis train requires that analysis libraries are compiled with a local AliRoot"); 
+      return kFALSE;
+   }   
+   Bool_t success = kTRUE;
+   // ROOT libraries
+   gSystem->Load("libTree.so");
+   gSystem->Load("libGeom.so");
+   gSystem->Load("libVMC.so");
+   gSystem->Load("libPhysics.so");
+   
+   // Load framework classes. Par option ignored here.
+   switch (imode) {
+      case 0:
+      case 2:
+            success &= LoadLibrary("libSTEERBase.so", mode);
+            success &= LoadLibrary("libESD.so", mode);
+            success &= LoadLibrary("libAOD.so", mode);
+            success &= LoadLibrary("libANALYSIS.so", mode);
+            success &= LoadLibrary("libANALYSISalice.so", mode);
+            success &= LoadLibrary("libCORRFW.so", mode);
+            gROOT->ProcessLine(".include $ALICE_ROOT/include");
+         break;
+      case 1:
+         Int_t ires = -1;
+         if (!gSystem->AccessPathName(AFversion)) ires = gProof->UploadPackage(AFversion);
+         if (ires < 0) {
+            success &= LoadLibrary("STEERBase", mode);
+            success &= LoadLibrary("ESD", mode);
+            success &= LoadLibrary("AOD", mode);
+            success &= LoadLibrary("ANALYSIS", mode);
+            success &= LoadLibrary("ANALYSISalice", mode);
+            success &= LoadLibrary("CORRFW", mode);
+         } else { 
+            ires = gProof->EnablePackage(AFversion);
+            if (ires<0) success = kFALSE;
+         }
+         break;         
+      default:
+         ::Error("LoadCommonLibraries", "Unknown run mode: %s", mode);
+         return kFALSE;
+   }
+   if (success) {
+      ::Info("LoadCommodLibraries", "Load common libraries:    SUCCESS");
+      ::Info("LoadCommodLibraries", "Include path for Aclic compilation:\n%s",
+              gSystem->GetIncludePath());
+   } else {           
+      ::Info("LoadCommodLibraries", "Load common libraries:    FAILED");
+   }   
+      
+   return success;
+}
+
+//______________________________________________________________________________
+TChain *CreateChain(const char *mode, const char *plugin_mode)
+{
+// Create the input chain
+   Int_t imode = -1;
+   if (!strcmp(mode, "LOCAL")) imode = 0;
+   if (!strcmp(mode, "PROOF")) imode = 1;
+   if (!strcmp(mode, "GRID"))  imode = 2;
+   TChain *chain = NULL;
+   // Local chain
+   switch (imode) {
+      case 0:
+            if (!strlen(dataset)) {
+               // Local ESD
+               chain = new TChain("esdTree");
+               if (gSystem->AccessPathName("/n60raid3/alice/roberta/MCBranch/AliESDs.root")) 
+                  ::Error("CreateChain", "File: AliESDs.root not in ./data dir");
+               else chain->Add("/n60raid3/alice/roberta/MCBranch/AliESDs.root");
+            } else {
+               // Interactive ESD
+               chain = CreateChainSingle(dataset, "esdTree");
+            }   
+         break;
+      case 1:
+         break;
+      case 2:
+            TString  treeName = "esdTree";
+            chain = CreateChainSingle("wn.xml", treeName);
+         break;      
+      default:   
+   }
+   if (chain && chain->GetNtrees()) return chain;
+   return NULL;
+}   
+
+//______________________________________________________________________________
+Bool_t LoadLibrary(const char *module, const char *mode, Bool_t rec=kFALSE)
+{
+// Load a module library in a given mode. Reports success.
+   Int_t imode = -1;
+   Int_t result;
+   TString smodule(module);
+   if (!strcmp(mode, "LOCAL")) imode = 0;
+   if (!strcmp(mode, "PROOF")) imode = 1;
+   if (!strcmp(mode, "GRID"))  imode = 2;
+   TString mod(module);
+   if (!mod.Length()) {
+      ::Error("LoadLibrary", "Empty module name");
+      return kFALSE;
+   }   
+   // If a library is specified, just load it
+   if (smodule.EndsWith(".so")) {
+      mod.Remove(mod.Index(".so"));
+      result = gSystem->Load(mod);
+      if (result < 0) {
+         ::Error("LoadLibrary", "Could not load library %s", module);
+         return kFALSE;
+      }
+      if (rec) anaLibs += Form("%s.so ",mod.Data()); 
+      return kTRUE;
+   } 
+   // Check if the library is already loaded
+   if (strlen(gSystem->GetLibraries(Form("%s.so", module), "", kFALSE)) > 0)
+      return kTRUE;    
+   switch (imode) {
+      case 0:
+      case 2:
+            result = gSystem->Load(Form("lib%s.so", module));
+            if (rec) anaLibs += Form("lib%s.so ", module);
+         break;
+      case 1:
+         result = gProof->UploadPackage(module);
+         if (result<0) {
+            result = gProof->UploadPackage(gSystem->ExpandPathName(Form("$ALICE_ROOT/%s.par", module)));
+            if (result<0) {
+               ::Error("LoadLibrary", "Could not find module %s.par in current directory nor in $ALICE_ROOT", module);
+               return kFALSE;
+            }
+         }   
+         result = gProof->EnablePackage(module);
+         break;
+      default:
+         return kFALSE;
+   }         
+   if (result < 0) {
+      ::Error("LoadLibrary", "Could not load module %s", module);
+      return kFALSE;
+   }
+   return kTRUE;
+}           
+
+//______________________________________________________________________________
+TChain* CreateChainSingle(const char* xmlfile, const char *treeName)
+{
+   printf("*******************************\n");
+   printf("*** Getting the ESD Chain   ***\n");
+   printf("*******************************\n");
+   TAlienCollection * myCollection  = TAlienCollection::Open(xmlfile);
+
+   if (!myCollection) {
+      ::Error("CreateChainSingle", "Cannot create an AliEn collection from %s", xmlfile) ;
+      return NULL ;
+   }
+
+   TChain* chain = new TChain(treeName);
+   myCollection->Reset() ;
+   while ( myCollection->Next() ) chain->Add(myCollection->GetTURL("")) ;
+   chain->ls();
+   return chain;
+}
+
+//______________________________________________________________________________
+Bool_t LoadAnalysisLibraries(const char *mode)
+{
+// Load common analysis libraries.
+   Bool_t success = kTRUE;
+      if (!LoadLibrary("PWG3base", mode, kTRUE) ||
+          !LoadLibrary("PWG3muon", mode, kTRUE)) return kFALSE;
+   ::Info("LoadAnalysisLibraries", "Load other libraries:   SUCCESS");
+   return kTRUE;
+}
+
diff --git a/PWG3/muon/CreateAlienHandler_FromESDToAOD.C b/PWG3/muon/CreateAlienHandler_FromESDToAOD.C
new file mode 100644 (file)
index 0000000..fb82c1f
--- /dev/null
@@ -0,0 +1,113 @@
+AliAnalysisGrid* CreateAlienHandler_FromESDToAOD()
+{
+//========================================================================
+// Macro to configure the GRID plugin
+// (see Alice offline web pages for definitions)
+//========================================================================
+
+// Check if user has a valid token, otherwise make one. This has limitations.
+// One can always follow the standard procedure of calling alien-token-init then
+//   source /tmp/gclient_env_$UID in the current shell
+//=====================================================================
+
+   if (!AliAnalysisGrid::CreateToken()) return NULL;
+   AliAnalysisAlien *plugin = new AliAnalysisAlien();
+
+// Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+//=====================================================================
+   plugin->SetRunMode("full");  // VERY IMPORTANT 
+
+// Set versions of used packages
+//=====================================================================
+   plugin->SetAPIVersion("V2.4");
+   plugin->SetROOTVersion("v5-23-04");
+   plugin->SetAliROOTVersion("v4-17-03");
+
+// Declare input data to be processed.
+//=====================================================================
+
+// Method 1: Create automatically XML collections using alien 'find' command.
+//===========
+// Define production directory LFN
+//   plugin->SetGridDataDir("/alice/cern.ch/user/a/arnaldi/FromESDToAOD/JPSI_generation/1001");
+// Set data search pattern
+//   plugin->SetDataPattern("*tag.root");
+// ...then add run numbers to be considered
+//   plugin->AddRunNumber(300000);
+//   plugin->AddRunNumber(1001);
+
+// Method 2: Declare existing data files (raw collections, xml collections, root file)
+//===========
+// If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir())
+// XML collections added via this method can be combined with the first method if
+// the content is compatible (using or not tags)
+   plugin->AddDataFile("/alice/cern.ch/user/a/arnaldi/FromESDToAOD/Plugin/FromESDToAOD.xml");
+
+// Define alien work directory where all files will be copied. Relative to alien $HOME.
+//=====================================================================
+   plugin->SetGridWorkingDir("FromESDToAOD/Plugin");
+   
+// Declare alien output directory. Relative to working directory.
+//=====================================================================
+   plugin->SetGridOutputDir("outputPlugin"); // In this case will be $HOME/work/output
+   
+// Declare the analysis source files names separated by blancs. 
+// Declare all libraries (other than the default ones for the framework). These will be
+// loaded by the generated analysis macro and compiled runtime.
+// Add par files, if needed.
+// Add all extra files (task .cxx/.h/.C) here.
+// (AddTaskTagCreation.C can be removed from SetAdditionalLibs, if available in the grid aliroot version)
+//=====================================================================
+   plugin->SetAdditionalLibs("libPWG3base.so libPWG3muon.so AddTaskTagCreation.C");
+
+// Declare the output file names separated by blancs.
+//=====================================================================
+// (can be like: file.root or file.root@ALICE::Niham::File)
+   plugin->SetOutputFiles("AliAODs.root AOD.tag.root");
+
+// Optionally define the files to be archived.
+//=====================================================================
+//   plugin->SetOutputArchive("log_archive.zip:stdout,stderr@ALICE::NIHAM::File root_archive.zip:*.root@ALICE::NIHAM::File");
+   plugin->SetOutputArchive("log_archive.zip:stdout,stderr");
+
+// Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+//=====================================================================
+   plugin->SetAnalysisMacro("analysisFromESDToAOD_Plugin.C");
+
+// Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+//=====================================================================
+   plugin->SetSplitMaxInputFileNumber(0);
+
+// Optionally set number of failed jobs that will trigger killing waiting sub-jobs.
+//=====================================================================
+   plugin->SetMaxInitFailed(5);
+
+// Optionally resubmit threshold.
+//=====================================================================
+   plugin->SetMasterResubmitThreshold(90);
+
+// Optionally set time to live (default 30000 sec)
+//=====================================================================
+   plugin->SetTTL(20000);
+
+// Optionally set input format (default xml-single)
+//=====================================================================
+   plugin->SetInputFormat("xml-single");
+
+// Optionally modify the name of the generated JDL (default analysis.jdl)
+//=====================================================================
+   plugin->SetJDLName("analysisFromESDToAOD_Plugin.jdl");
+
+// Optionally modify job price (default 1)
+//=====================================================================
+   plugin->SetPrice(1); 
+
+// Optionally modify split mode (default 'se')    
+//=====================================================================
+   plugin->SetSplitMode("se");
+
+// Optionally define preferred SE
+//=====================================================================
+   plugin->SetPreferedSE("ALICE::Torino::DPM");
+   return plugin;
+} 
diff --git a/PWG3/muon/ReadAOD_MCBranch.C b/PWG3/muon/ReadAOD_MCBranch.C
new file mode 100644 (file)
index 0000000..1829a9d
--- /dev/null
@@ -0,0 +1,179 @@
+void ReadAOD_MCBranch(char *filein="AliAODs.root"){
+
+//=============================================================================
+// Macro to read AOD+MC Branch: 
+// in this example the macro reads generated and reconstructed JPsi 
+// Plots refer to all MC generated particles, to MC particles which have
+// been reconstructed and reconstructed tracks
+//
+// Each aod event has 3 MC particles (J/Psi, mu1, mu2)
+// To access MC muons:   mctrackall->IsPhysicalPrimary()
+// To access JPsi muons: mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()
+
+// AliStack::IsPhysicalPrimary() 
+// Test if a particle is a physical primary according to the following definition:
+// Particles produced in the collision including products of strong and
+// electromagnetic decay and excluding feed-down from weak decays of strange
+// particles.
+//=============================================================================
+
+  gStyle->SetFrameFillColor(0);
+  gStyle->SetCanvasColor(0);
+  gStyle->SetPalette(1); 
+
+// Open AOD file    
+  TFile *file = new TFile(filein);
+  TTree *aodTree = file->Get("aodTree");
+  
+// Read AOD event    
+  AliAODEvent *event = new AliAODEvent();
+  event->ReadFromTree(aodTree);
+  
+// Read MC info   
+  TClonesArray *mcarray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
+  if(!mcarray) {
+    printf("No MC info in AOD!\n");
+    return;
+  }  
+  
+// Book Histos
+  TH2D *hpt_mc_rec = new TH2D("hpt_mc_rec","Pt correlations MC-REC",50,0.,10.,50,0.,10.);
+  TH1D *hpt_rec = new TH1D("hpt_rec","Pt REC tracks",20,0.,10.);
+  TH1D *hpt_mc = new TH1D("hpt_mc","PT MC particles corresponding to rec tracks",20,0.,10.);
+  TH1D *hpt_allmc = new TH1D("hpt_allmc","PT all MC particles",20,0.,10.);
+  TH1D *hpt_primarymc = new TH1D("hpt_primarymc","PT MC primaries (J/Psi and muons)",20,0.,10.);
+  TH1D *hpt_physprimarymc = new TH1D("hpt_physprimarymc","PT Muons MC",20,0.,10.);
+  TH1D *hpt_notphysprimarymc = new TH1D("hpt_notphysprimarymc","Pt JPsi MC",20,0.,10.);
+  TH1D *hpt_diff = new TH1D("hpt_diff","Pt: Rec - MC",100,-5.,5.);
+
+  TH1D *hCharge_primarymc = new TH1D("hCharge_primarymc","Charge MC primary particle",10,-5.,5.);
+  TH1D *hCharge_physprimmc = new TH1D("hCharge_physprimmc","Charge Muons MC",10,-5.,5.);
+  TH1D *hCharge_notphysprimmc = new TH1D("hCharge_notphysprimmc","Charge JPsi MC",10,-5.,5.);
+  TH1D *hCharge_rec = new TH1D("hCharge_rec","hCharge_rec",10,-5.,5.);
+  TH1D *hCharge_mc = new TH1D("hCharge_mc","hCharge_mc",10,-5.,5.);
+
+  TH1D *hMass_mc = new TH1D("hMass_mc","Mass MC particles",100,0.,5..);
+  TH1D *hMass_rec = new TH1D("hMass_rec","Mass rec tracks",100,0.,5..);
+  TH1D *hMass_allmc = new TH1D("hMass_allmc","Mass all MC particles",100,0.,5..);
+  TH1D *hMass_primarymc = new TH1D("hMass_primarymc","Mass primary particles MC",100,0.,5..);
+  TH1D *hMass_physprimarymc = new TH1D("hMass_physprimarymc","Mass muons MC",100,0.,5..);
+  TH1D *hMass_notphysprimarymc = new TH1D("hMass_notphysprimarymc","Mass JPsi MC",100,0.,5..);
+  TH1D *hMassJPsi_mc = new TH1D("hMassJPsi_mc","Mass JPsi MC",100,0.,5..);
+  
+// Loop on AOD events
+  Int_t nev=0;
+
+  while(aodTree->GetEvent(nev++)){  
+
+// Loop on all MC tracks    
+    for(int ii=0;ii<mcarray->GetEntries();ii++){
+      AliAODMCParticle *mctrackall = (AliAODMCParticle*) mcarray->At(ii);     
+      // all particles
+      hpt_allmc->Fill(mctrackall->Pt());
+      hMass_allmc->Fill(mctrackall->M());
+      
+      // Muons
+      if(mctrackall->IsPhysicalPrimary()) {
+        hCharge_physprimmc->Fill(mctrackall->Charge());
+        hpt_physprimarymc->Fill(mctrackall->Pt());
+        hMass_physprimarymc->Fill(mctrackall->M());
+      }
+       
+      if(mctrackall->IsPrimary()) {
+        hpt_primarymc->Fill(mctrackall->Pt());
+        hMass_primarymc->Fill(mctrackall->M());
+       hCharge_primarymc->Fill(mctrackall->Charge());
+      }
+        
+      // JPsi
+      if(mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()) {
+        hCharge_notphysprimmc->Fill(mctrackall->Charge());
+        hMass_notphysprimarymc->Fill(mctrackall->M());
+        hpt_notphysprimarymc->Fill(mctrackall->Pt());
+      }        
+    }
+    
+// Loop on reconstructed tracks    
+    for(int i=0;i<event->GetNumberOfTracks();i++){
+      AliAODTrack *track = dynamic_cast<AliAODTrack*>(event->GetTrack(i));
+      if(!track) {
+        printf("Error: no track! \n");
+        continue;
+      } 
+      hpt_rec->Fill(track->Pt());
+      hMass_rec->Fill(track->M());
+      hCharge_rec->Fill(track->Charge());
+      
+// Read MC particle corresponding to reconstructed tracks   
+      AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(track->GetLabel());     
+      if(!mctrack){
+        printf("Error: no MC track! \n");
+       continue;
+      }
+      //printf("nev=%d i=%d mc=%f rec=%f\n",nev-1,i,mctrack->Pt(),track->Pt());
+      hpt_mc->Fill(mctrack->Pt());
+      hpt_mc_rec->Fill(track->Pt(),mctrack->Pt());
+      hpt_diff->Fill(track->Pt()-mctrack->Pt());
+      hMass_mc->Fill(mctrack->M());
+      hCharge_mc->Fill(mctrack->Charge());
+   }
+ }
+// Plots 
+ TCanvas *c = new TCanvas("c","PT spectra",20,20,600,600);
+ c->Divide(3,3);
+ c->cd(1);
+ hpt_primarymc->Draw();
+ hpt_primarymc->GetXaxis()->SetTitle("Pt MC primary particles");
+ c->cd(2);  
+ hpt_mc->GetXaxis()->SetTitle("Pt MC particles");
+ hpt_mc->Draw();
+ hpt_mc->SetLineColor(4);
+ c->cd(3);
+ hpt_rec->GetXaxis()->SetTitle("Rec Pt track");
+ hpt_rec->Draw();
+ hpt_rec->SetLineColor(2);
+ c->cd(4);
+ hpt_physprimarymc->Draw();
+ hpt_physprimarymc->GetXaxis()->SetTitle("Pt MC muon track");
+ c->cd(5);
+ hpt_notphysprimarymc->Draw();
+ hpt_notphysprimarymc->GetXaxis()->SetTitle("Pt MC JPsi track");
+ c->cd(7);
+ hpt_mc_rec->Draw("colz");
+ hpt_mc_rec->GetXaxis()->SetTitle("Pt rec track");
+ hpt_mc_rec->GetYaxis()->SetTitle("Pt MC track");
+ c->cd(8);
+ hpt_diff->Draw();
+ hpt_diff->GetXaxis()->SetTitle("Rec-MC Pt");
+
+ TCanvas *c1 = new TCanvas("c1","Mass distributions",40,40,600,600);
+ c1->Divide(2,2);
+ c1->cd(1);
+ hMass_primarymc->Draw();
+ hMass_primarymc->GetXaxis()->SetTitle("Mass MC primary particles");
+ c1->cd(2);
+ hMass_physprimarymc->Draw();
+ hMass_physprimarymc->GetXaxis()->SetTitle("Mass MC muon particles");
+ c1->cd(3);
+ hMass_notphysprimarymc->Draw();
+ hMass_notphysprimarymc->GetXaxis()->SetTitle("Mass MC J/Psi");
+ c1->cd(4);
+ hMass_rec->Draw();
+ hMass_rec->GetXaxis()->SetTitle("Mass Rec tracks");
+
+ TCanvas *c2= new TCanvas("c2","Charge distributions",60,60,600,600);
+ c2->Divide(2,2);
+ c2->cd(1);
+ hCharge_primarymc->Draw();
+ hCharge_primarymc->GetXaxis()->SetTitle("Charge MC primary particles");
+ c2->cd(2);
+ hCharge_physprimmc->Draw();
+ hCharge_physprimmc->GetXaxis()->SetTitle("Charge muon MC");;
+ c2->cd(3);
+ hCharge_notphysprimmc->Draw();
+ hCharge_notphysprimmc->GetXaxis()->SetTitle("Charge J/Psi MC");;
+ c2->cd(4);
+ hCharge_rec->Draw();
+ hCharge_rec->GetXaxis()->SetTitle("Charge rec tracks");
+}