PWG4/anaPhos.C

   1 Bool_t gIsAnalysisLoaded = kFALSE ;
   2
   3 //______________________________________________________________________
   4 Bool_t LoadLib( const char* pararchivename)
   5 {
   6   // Loads the AliRoot required libraries from a tar file
   7
   8   Bool_t rv = kTRUE ;
   9
  10   char cdir[1024] ;
  11   sprintf(cdir, "%s", gSystem->WorkingDirectory() ) ;
  12
  13   // Setup par File
  14   if (pararchivename) {
  15    char parpar[80] ;
  16    sprintf(parpar, "%s.par", pararchivename) ;
  17    if ( gSystem->AccessPathName(parpar) ) {
  18     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
  19     char processline[1024];
  20     sprintf(processline, ".! make %s", parpar) ;
  21     cout << processline << endl ;
  22     gROOT->ProcessLine(processline) ;
  23     gSystem->ChangeDirectory(cdir) ;
  24     sprintf(processline, ".! mv /tmp/%s .", parpar) ;
  25     gROOT->ProcessLine(processline) ;
  26     sprintf(processline,".! tar xvzf %s",parpar);
  27     gROOT->ProcessLine(processline);
  28    }
  29    gSystem->ChangeDirectory(pararchivename);
  30
  31     // check for BUILD.sh and execute
  32     if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
  33       printf("*** Building PAR archive  %s  ***\n", pararchivename);
  34
  35       if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
  36         printf("Cannot Build the PAR Archive %s! - Abort!", pararchivename) );
  37         return kFALSE ;
  38       }
  39     }
  40
  41     // check for SETUP.C and execute
  42     if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
  43       printf("*** Setup PAR archive  %s     ***\n", pararchivename);
  44       gROOT->Macro("PROOF-INF/SETUP.C");
  45     }
  46   }
  47
  48
  49   if ( strstr(pararchivename, "AnalysisCheck") ) {
  50     gSystem->Load("libSpectrum.so");
  51   }
  52
  53   printf("lib%s done\n", pararchivename);
  54
  55   gSystem->ChangeDirectory(cdir);
  56
  57   gIsAnalysisLoaded = kTRUE ;
  58   return rv ;
  59 }
  60
  61 //______________________________________________________________________
  62 void anaPhos(const Int_t kEvent=10)
  63 {
  64   if (! gIsAnalysisLoaded ) {
  65     LoadLib("ESD") ;
  66     LoadLib("AOD") ;
  67     LoadLib("ANALYSIS") ;
  68     LoadLib("PWG4Gamma") ;
  69   }
  70
  71   // create the analysis goodies object
  72   AliAnalysisGoodies * ag = new AliAnalysisGoodies() ;
  73
  74   // definition of analysis tasks
  75
  76   // first task
  77   AliAnaGammaPhos * phostask = new AliAnaGammaPhos("GammaPhos") ;
  78   ag->ConnectInput(phostask, TChain::Class(), 0) ;
  79   ag->ConnectOuput(phostask, TTree::Class(), 0, "AOD") ;
  80   AliAnalysisDataContainer * outGammaPhos = ag->ConnectOuput(phostask, TList::Class(), 1) ;
  81
  82   AliAnaScale * scale = new AliAnaScale("ScaledGammaPhos") ;
  83   ag->ConnectInput(scale, outGammaPhos, 0) ;
  84   ag->ConnectOuput(scale, TList::Class(), 0) ;
  85
  86
  87   // get the data to analyze
  88
  89   // definition of Tag cuts
  90   const char * runCuts = 0x0 ;
  91   const char * evtCuts = 0x0 ;
  92   const char * lhcCuts = 0x0 ;
  93   const char * detCuts = 0x0 ;
  94
  95 //"fEventTag.fNPHOSClustersMin == 1 && fEventTag.fNEMCALClustersMin == 1" ;
  96
  97
  98   TString input = gSystem->Getenv("ANA_INPUT") ;
  99   if ( input != "") {
 100     char argument[1024] ;
 101     if ( input.Contains("tag?") ) {
 102       //create the ESD collection from the tag collection
 103       input.ReplaceAll("tag?", "") ;
 104       const char * collESD = "esdCollection.xml" ;
 105       ag->MakeEsdCollectionFromTagCollection(runCuts, lhcCuts, detCuts, evtCuts, input.Data(), collESD) ;
 106       sprintf(argument, "esd?%s", collESD) ;
 107     } else if ( input.Contains("esd?") )
 108       sprintf(argument, "%s", input.Data()) ;
 109     ag->Process(argument) ;
 110
 111   } else {
 112
 113     TChain* analysisChain = new TChain("esdTree") ;
 114     //   input = "alien:///alice/cern.ch/user/a/aliprod/prod2006_2/output_pp/105/411/AliESDs.root" ;
 115     //   analysisChain->AddFile(input);
 116     input = "AliESDs.root" ;
 117     const char * kInDir = gSystem->Getenv("OUTDIR") ;
 118     if ( kInDir ) {
 119       if ( ! gSystem->cd(kInDir) ) {
 120         printf("%s does not exist\n", kInDir) ;
 121         return ;
 122       }
 123       Int_t event, skipped=0 ;
 124       char file[120] ;
 125       Double_t xsection = 0., ntrials = 0. ;
 126       for (event = 0 ; event < kEvent ; event++) {
 127         sprintf(file, "%s/%d/AliESDs.root", kInDir,event) ;
 128         TFile * fESD = 0 ;
 129         if ( fESD = TFile::Open(file))
 130           if ( fESD->Get("esdTree") ) {
 131             printf("++++ Adding %s\n", file) ;
 132             analysisChain->AddFile(file);
 133             Double_t * rv = ReadXsection(kInDir, event) ;
 134             xsection = xsection + rv[0] ;
 135             ntrials  = ntrials  + rv[1] ;
 136             cout << xsection << " " << ntrials << endl ;
 137
 138           }
 139           else {
 140             printf("---- Skipping %s\n", file) ;
 141             skipped++ ;
 142           }
 143       }
 144       printf("number of entries # %lld, skipped %d\n", analysisChain->GetEntries(), skipped*100) ;
 145     }
 146     else
 147       analysisChain->AddFile(input);
 148
 149     scale->Set(xsection/ntrials) ;
 150     ag->Process(analysisChain) ;
 151   }
 152   return ;
 153 }
 154
 155 //______________________________________________________________________
 156 void Merge(const char * xml, const char * sub, const char * out)
 157 {
 158   if (! gIsAnalysisLoaded )
 159     LoadLib("ESD") ;
 160
 161   AliAnalysisGoodies * ag = new AliAnalysisGoodies() ;
 162   ag->Merge(xml, sub, out) ;
 163 }
 164
 165 Double_t * ReadXsection(const char * inDir, const Int_t event)
 166 {
 167   // Read the PYTHIA statistics from the file pyxsec.root created by
 168   // the function WriteXsection():
 169   // integrated cross section (xsection) and
 170   // the  number of Pyevent() calls (ntrials)
 171   // and calculate the weight per one event xsection/ntrials
 172   // The spectrum calculated by a user should be
 173   // multiplied by this weight, something like this:
 174   // TH1F *userSpectrum ... // book and fill the spectrum
 175   // userSpectrum->Scale(weight)
 176   //
 177   // Yuri Kharlov 19 June 2007
 178
 179   Double_t xsection;
 180   UInt_t    ntrials;
 181
 182   char cfile[80] ;
 183   sprintf(cfile, "%s/%d/pyxsec.root", inDir, event) ;
 184
 185   TFile *file = new TFile(cfile,"readonly");
 186   if ( ! file ) {
 187     AliFatal(Form("could not open %s", cfile)) ;
 188     exit(1) ;
 189   }
 190   TTree *tree = file->Get("Xsection");
 191   tree->SetBranchAddress("xsection",&xsection);
 192   tree->SetBranchAddress("ntrials",&ntrials);
 193   tree->GetEntry(0);
 194   cout << "Cross section = "<<xsection<<" mb, N trials = "<<ntrials<<endl;
 195   Double_t rv[2] ;
 196   rv[0] = xsection ;
 197   rv[1] = ntrials ;
 198   return rv ;
 199 }