automatic histogram scaling for AODs now available
authorkread <kread@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Sep 2009 06:13:46 +0000 (06:13 +0000)
committerkread <kread@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Sep 2009 06:13:46 +0000 (06:13 +0000)
PWG4/macros/electrons/ConfigJetAnalysisFastJet.C
PWG4/macros/electrons/anaJete.C
PWG4/macros/electrons/anaJete.sh

index 8ab5dd4..11d3d99 100755 (executable)
@@ -17,9 +17,9 @@ AliJetFinder*  ConfigJetAnalysis()
     printf("========================== \n");\r
 \r
     Bool_t kInputIsESD = kTRUE;     //uncomment for input ESD\r
-    //Bool_t kInputIsESD = kFALSE;    //uncomment for input AODs\r
+  //Bool_t kInputIsESD = kFALSE;    //uncomment for input AODs\r
     Bool_t kFollowsFilter = kTRUE;  //uncomment if follows ESD filter task\r
-    //Bool_t kFollowsFilter = kFALSE; //uncomment if no ESD filter task\r
+  //Bool_t kFollowsFilter = kFALSE; //uncomment if no ESD filter task\r
 \r
 \r
     // Define the grids\r
index 03dad04..62c116a 100755 (executable)
@@ -53,6 +53,10 @@ char sconfig1[1024] = "ConfigPWG4AODtoAOD";        //"ConfigAnalysis";
 char sconfig2[1024] = "ConfigJetAnalysisFastJet.C";//"ConfigAnalysis";\r
 char sconfig3[1024] = "ConfigAnalysisElectron";    //"ConfigAnalysis";\r
 \r
+//Initialize the cross section and ntrials values. Do not modify.\r
+Double_t xsection = 0;\r
+Float_t ntrials = 0;\r
+\r
 void anaJete()\r
 {\r
   // Main\r
@@ -84,6 +88,7 @@ void anaJete()
 \r
   if (mode==mLocal || mode==mLocalCAF || mode == mGRID) {\r
     CreateChain(mode, chain, chainxs);  \r
+    cout<<"Chain created"<<endl;\r
   }\r
 \r
   if( chain || mode==mPROOF || mode==mPLUGIN ){\r
@@ -167,20 +172,22 @@ void anaJete()
       AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Loose");\r
       //esdTrackCutsL->SetMinNClustersTPC(50);\r
       //esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5);\r
-      //esdTrackCutsL->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);\r
       //esdTrackCutsL->SetRequireTPCRefit(kTRUE);\r
-      //esdTrackCutsL->SetMinNsigmaToVertex(3); //keep commented out\r
-      //esdTrackCutsL->SetRequireSigmaToVertex(kTRUE); //keep commented out\r
+      //esdTrackCutsL->SetMaxDCAToVertexXY(2.4);\r
+      //esdTrackCutsL->SetMaxDCAToVertexZ(3.2);\r
+      //esdTrackCutsL->SetDCAToVertex2D(kTRUE);\r
+      //esdTrackCutsL->SetRequireSigmaToVertex(kFALSE);\r
       //esdTrackCutsL->SetAcceptKinkDaughters(kFALSE);\r
       //\r
       // hard\r
       AliESDtrackCuts* esdTrackCutsH = new AliESDtrackCuts("AliESDtrackCuts", "Hard");\r
       esdTrackCutsH->SetMinNClustersTPC(100);\r
       esdTrackCutsH->SetMaxChi2PerClusterTPC(2.0);\r
-      esdTrackCutsH->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);\r
       esdTrackCutsH->SetRequireTPCRefit(kTRUE);\r
-      //esdTrackCutsH->SetMinNsigmaToVertex(2);\r
-      //esdTrackCutsH->SetRequireSigmaToVertex(kTRUE);\r
+      esdTrackCutsH->SetMaxDCAToVertexXY(2.4);\r
+      esdTrackCutsH->SetMaxDCAToVertexZ(3.2);\r
+      esdTrackCutsH->SetDCAToVertex2D(kTRUE);\r
+      esdTrackCutsH->SetRequireSigmaToVertex(kFALSE);\r
       esdTrackCutsH->SetAcceptKinkDaughters(kFALSE);\r
       //\r
       AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");\r
@@ -239,30 +246,31 @@ void anaJete()
     //------------------------  \r
     //Scaling task\r
     //-----------------------\r
-    Int_t nfiles = chainxs->GetEntries();\r
+    cout<<">>> Scaling Task"<<endl;\r
+    Int_t nfiles = chain->GetListOfFiles()->GetEntriesFast();//chainxs->GetEntries();\r
     Int_t nevents = chain->GetEntries();\r
     cout<<"Get? "<<kGetXSectionFromFileAndScale<<" nfiles "<<nfiles<<" nevents "<<nevents<<endl;\r
     if(kGetXSectionFromFileAndScale && nfiles > 0){\r
-      //cout<<"Init AnaScale"<<endl;\r
-      //Get the cross section\r
-      Double_t xsection=0; \r
-      Float_t ntrials = 0;\r
-      GetAverageXsection(chainxs, xsection, ntrials);\r
-      \r
+      cout<<"Init AnaScale"<<endl;\r
       AliAnaScale * scale = new AliAnaScale("scale") ;\r
-\r
-      cout<<"Scale factor "<<xsection/(ntrials/kNumberOfEventsPerFile)/nevents<<"  "<<xsection/ntrials/nfiles<<endl;\r
-      scale->Set(xsection/ntrials/nfiles) ;\r
+      cout<<"Summed xsection "<<xsection<<" Summed ntrials "<<ntrials<<" total events "<<nevents<<endl;\r
+      //Calculate the average\r
+      xsection/=nfiles;\r
+      ntrials/=nfiles;\r
+      cout<<"Average xsection "<<xsection<<" Average ntrials "<<ntrials<<" total events "<<nevents<<endl;\r
+      Double_t scaleFactor =   xsection/ntrials/nevents ;\r
+      cout<<"Scale factor "<<scaleFactor<<endl;\r
+      scale->Set(scaleFactor) ;\r
       scale->MakeSumw2(kTRUE);//If you want histograms with error bars set to kTRUE\r
       scale->SetDebugLevel(2);\r
       mgr->AddTask(scale);\r
       \r
-      AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("histosscaled", TList::Class(),\r
-                                                               AliAnalysisManager::kOutputContainer, "histosscaled.root");\r
+      AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("histosscaled", \r
+                  TList::Class(), AliAnalysisManager::kOutputContainer, "histosscaled.root");\r
       mgr->ConnectInput  (scale,     0, coutput2);\r
       mgr->ConnectOutput (scale,     0, coutput3 );\r
     }\r
-    \r
+   \r
     //-----------------------\r
     // Run the analysis\r
     //-----------------------    \r
@@ -284,7 +292,7 @@ void anaJete()
     else if (mode==mPLUGIN)\r
       mgr->StartAnalysis(smode.Data());\r
     else\r
-      mgr->StartAnalysis(smode.Data(),chain,20);\r
+      mgr->StartAnalysis(smode.Data(),chain);\r
 \r
     cout <<" Analysis ended sucessfully "<< endl ;\r
 \r
@@ -451,6 +459,24 @@ void SetupPar(char* pararchivename)
 \r
 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){\r
   //Fills chain with data\r
+\r
+  TString datafile="";\r
+  if(kInputData == "ESD") datafile = "AliESDs.root" ;\r
+  else if(kInputData == "AOD") datafile = "AliAOD.root" ;\r
+  else if(kInputData == "MC")  datafile = "galice.root" ;\r
+\r
+  if(kInputData == "AOD") kXSFileName = "pyxsec_hists.root";\r
+\r
+  char * kGener = gSystem->Getenv("GENER");\r
+  if(kGener) {\r
+    cout<<"GENER "<<kGener<<endl;\r
+    if(!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";\r
+    else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";\r
+    else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;\r
+  }\r
+  else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;\r
+\r
+\r
   TString ocwd = gSystem->WorkingDirectory();\r
   \r
   //-----------------------------------------------------------\r
@@ -500,26 +526,13 @@ void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
       //if(gSystem->Getenv("XSFILE"))  \r
       //kXSFileName = gSystem->Getenv("XSFILE") ; \r
       //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;  \r
-      char * kGener = gSystem->Getenv("GENER");\r
-      if(kGener) {\r
-       cout<<"GENER "<<kGener<<endl;\r
-       if(!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";\r
-       else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";\r
-       else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;\r
-      }\r
-\r
-      else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;\r
 \r
       cout<<"INDIR : "<<kInDir<<endl;\r
       cout<<"NFILES : "<<kFile<<endl;\r
       cout<<"PATTERN: " <<kPattern<<endl;\r
       cout<<"XSFILE  : "<<kXSFileName<<endl;\r
       \r
-      TString datafile="";\r
-      if(kInputData == "ESD") datafile = "AliESDs.root" ;\r
-      else if(kInputData == "AOD") datafile = "AliAOD.root" ;\r
-      else if(kInputData == "MC")  datafile = "galice.root" ;\r
-      \r
+\r
       //Loop on ESD files, add them to chain\r
       Int_t event =0;\r
       Int_t skipped=0 ; \r
@@ -529,13 +542,14 @@ void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
       for (event = 0 ; event < kFile ; event++) {\r
        sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; \r
        sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ; \r
-       TFile * fESD = 0 ; \r
+       TFile * dataFile = 0 ; \r
        //Check if file exists and add it, if not skip it\r
-       if ( fESD = TFile::Open(file)) {\r
-         if ( fESD->Get(kTreeName) ) { \r
-           printf("++++ Adding %s\n", file) ;\r
+       if ( dataFile = TFile::Open(file)) {\r
+         if ( dataFile->Get(kTreeName) ) { \r
+           Int_t nEventsPerFile = ((TTree*) dataFile->Get(kTreeName)) ->GetEntries();\r
+           printf("++++ Adding %s, with %d events \n", file, nEventsPerFile) ;\r
            chain->AddFile(file);\r
-           if(kGetXSectionFromFileAndScale)chainxs->Add(filexs) ; \r
+           if(kGetXSectionFromFileAndScale) GetXsection(nEventsPerFile, filexs);       \r
          }\r
        }\r
        else { \r
@@ -587,12 +601,23 @@ void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
    \r
     // Makes the ESD chain \r
     printf("*** Getting the Chain       ***\n");\r
+    Int_t nEventsPerFile = 0;\r
     for (Int_t index = 0; index < result->GetEntries(); index++) {\r
       TString alienURL = result->GetKey(index, "turl") ; \r
       cout << "================== " << alienURL << endl ; \r
       chain->Add(alienURL) ; \r
-      alienURL.ReplaceAll("AliESDs.root",kXSFileName);\r
-      chainxs->Add(alienURL) ; \r
+\r
+      if(kGetXSectionFromFileAndScale){\r
+       //Get the number of events per file.\r
+        //Do it only once, no need to open all the files.\r
+        if(i ==0 ) {\r
+          TFile * df = TFile::Open(alienURL);\r
+          nEventsPerFile = ((TTree*) df->Get(kTreeName)) ->GetEntries();\r
+          dataFile->Close();\r
+        } \r
+        alienURL.ReplaceAll(datafile,kXSFileName);\r
+        GetXsection(nEventsPerFile, alienURL);//chainxs->Add(alienURL) ; \r
+      }\r
     }\r
   }// xml analysis\r
 \r
@@ -600,64 +625,56 @@ void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
 }\r
 \r
 //________________________________________________\r
-void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)\r
+void GetXsection(Int_t nEventsPerFile, TString filexs)\r
 {\r
-  // Read the PYTHIA statistics from the file pyxsec.root created by\r
-  // the function WriteXsection():\r
-  // integrated cross section (xsection) and\r
-  // the  number of Pyevent() calls (ntrials)\r
-  // and calculate the weight per one event xsection/ntrials\r
-  // The spectrum calculated by a user should be\r
-  // multiplied by this weight, something like this:\r
-  // TH1F *userSpectrum ... // book and fill the spectrum\r
-  // userSpectrum->Scale(weight)\r
-  //\r
-  // Yuri Kharlov 19 June 2007\r
-  // Gustavo Conesa 15 April 2008\r
-  Double_t xsection = 0;\r
-  UInt_t    ntrials = 0;\r
-  xs = 0;\r
-  ntr = 0;\r
-  \r
-  Int_t nfiles =  tree->GetEntries()  ;\r
-  if (tree && nfiles > 0) {\r
-    tree->SetBranchAddress("xsection",&xsection);\r
-    tree->SetBranchAddress("ntrials",&ntrials);\r
-    for(Int_t i = 0; i < nfiles; i++){\r
-      tree->GetEntry(i);\r
-      xs += xsection ;\r
-      ntr += ntrials ;\r
-      cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl; \r
-    }\r
-    \r
-    xs =   xs /  nfiles;\r
-    ntr =  ntr / nfiles;\r
-    cout << "-----------------------------------------------------------------"<<endl;\r
-    cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl; \r
-    cout << "-----------------------------------------------------------------"<<endl;\r
-  } \r
-  else cout << " >>>> Empty tree !!!! <<<<< "<<endl;\r
-  \r
+  // Get the cros section from the corresponding file in the directory\r
+  // where the data sits.\r
+  // The xsection and ntrials global variables are updated per each file.\r
+  // The average of these cuantities should be calculated after.\r
+\r
+  TFile *fxs = new TFile(filexs,"R");\r
+  if(kInputData =="AOD") { //needs improvement, in case of train with ESDs, reading output AODs for example this is wrong.\r
+    TList *l = (TList*) fxs->Get("cFilterList");\r
+    TH1F * hxs = l->FindObject("h1Xsec") ;\r
+    TH1F * htrial = l->FindObject("h1Trials") ;\r
+    if(htrial->GetEntries()!=hxs->GetEntries() || htrial->GetEntries()==0){\r
+      cout<<"Careful!!! Entries in histo for cross section "<<hxs->GetEntries()<< ", for trials "<<htrial->GetEntries()<<endl;\r
+      continue;\r
+    }          \r
+    xsection += hxs->GetBinContent(1);\r
+    ntrials  += htrial->GetBinContent(1)/htrial->GetEntries()/nEventsPerFile;\r
+    cout << "Chain: xsection " <<hxs->GetBinContent(1)<<" ntrials "<<htrial->GetBinContent(1)/htrial->GetEntries()<<endl; \r
+  }\r
+  else {\r
+    Double_t xs = 0;\r
+    UInt_t ntr = 0;\r
+    TTree * xstree = (TTree*)fxs->Get("Xsection");\r
+    xstree->SetBranchAddress("xsection",&xs);\r
+    xstree->SetBranchAddress("ntrials",&ntr);\r
+    xstree->GetEntry(0);\r
+    cout << "Chain: xsection " <<xs<<" ntrials "<<ntr<<endl;\r
+    xsection += xs ;\r
+    ntrials += ntr/nEventsPerFile;     \r
+  }  \r
 }\r
 \r
 void ProcessEnvironment(){\r
 \r
-       if (gSystem->Getenv("MODE"))\r
-               mode = atoi(gSystem->Getenv("MODE"));\r
+  if (gSystem->Getenv("MODE"))\r
+     mode = atoi(gSystem->Getenv("MODE"));\r
 \r
-       if (gSystem->Getenv("CONFIG1"))\r
-               sprintf(sconfig1,gSystem->Getenv("CONFIG1"));\r
+  if (gSystem->Getenv("CONFIG1"))\r
+     sprintf(sconfig1,gSystem->Getenv("CONFIG1"));\r
 \r
-       if (gSystem->Getenv("CONFIG2"))\r
-               sprintf(sconfig2,gSystem->Getenv("CONFIG2"));\r
+  if (gSystem->Getenv("CONFIG2"))\r
+      sprintf(sconfig2,gSystem->Getenv("CONFIG2"));\r
 \r
-       if (gSystem->Getenv("CONFIG3"))\r
-               sprintf(sconfig3,gSystem->Getenv("CONFIG3"));\r
+  if (gSystem->Getenv("CONFIG3"))\r
+      sprintf(sconfig3,gSystem->Getenv("CONFIG3"));\r
 \r
-       if (gSystem->Getenv("SEVENT"))\r
-               sevent = atoi (gSystem->Getenv("SEVENT"));\r
+  if (gSystem->Getenv("SEVENT"))\r
+      sevent = atoi (gSystem->Getenv("SEVENT"));\r
        \r
-       printf("PROCESS: Variables: mode %d, config1 %s, config2 %s, config3 %s, sevent %d\n", mode,sconfig1,sconfig2,sconfig3,sevent);\r
+      printf("PROCESS: Variables: mode %d, config1 %s, config2 %s, config3 %s, sevent %d\n", mode,sconfig1,sconfig2,sconfig3,sevent);\r
 \r
-       \r
 }\r
index c561844..3b27d76 100755 (executable)
@@ -10,7 +10,6 @@ export PATTERN=
 #export PATTERN=01
 export NFILES=1
 export MODE=0
-export CONFIG1=ConfigPWG4AODtoAOD
 export CONFIG2=ConfigJetAnalysisFastJet.C
 export CONFIG3=ConfigAnalysisElectron
 export SIMPATH=/work2/data/bjetfilter/AOD/117005/004