C++ TObjArray::FindObject方法代码示例

void fitter(
            const string workDirName="Test", // Working directory
            // Select the type of datasets to fit
            bool fitData     = true,         // Fits Data if true, otherwise fits MC
            bool fitPbPb     = true,         // Fits PbPb datasets
            bool fitPP       = true,         // Fits PP datasets
            // Select the type of object to fit
            bool incJpsi     = true,         // Includes Jpsi model
            bool incPsi2S    = true,         // Includes Psi(2S) model
            bool incBkg      = true,         // Includes Background model
            // Select the fitting options
            bool cutCtau     = false,        // Apply prompt ctau cuts
            bool doSimulFit  = false,        // Do simultaneous fit
            bool wantPureSMC = false,        // Flag to indicate if we want to fit pure signal MC
            int  numCores    = 2,            // Number of cores used for fitting
            // Select the drawing options
            bool setLogScale = true,         // Draw plot with log scale
            bool incSS       = false,        // Include Same Sign data
            bool zoomPsi     = false,        // Zoom Psi(2S) peak on extra pad
            int  nBins       = 54            // Number of bins used for plotting
            ) 
{
  // -------------------------------------------------------------------------------
  // STEP 0: INITIALIZE THE FITTER WORK ENVIROMENT
  // The work enviroment is divided as follows:
  /*
    main |-> Macros: Contain all the macros
         |-> Input   |-> <WorkDir> : Contain Input File, Bin and Parameter List for a given work directory (e.g. 20160201)
	 |-> Output  |-> <WorkDir> : Contain Output Plots and Results for a given work directory (e.g. 20160201)
	 |-> DataSet : Contain all the datasets (MC and Data)
  */

  if (!checkSettings(fitData, fitPbPb, fitPP, incJpsi, incPsi2S, incBkg, cutCtau, doSimulFit, wantPureSMC, setLogScale, zoomPsi, incSS, numCores, nBins)) { return; }

  map<string,string> DIR;
  if(!iniWorkEnv(DIR, workDirName)){ return; }
 
  // -------------------------------------------------------------------------------
  // STEP 1: CREATE/LOAD THE ROODATASETS
  /*
    Input : List of TTrees with format:  TAG <tab> FILE_NAME
    Output: Collection of RooDataSets splitted by tag name, including OS and SS dimuons.
  */
  
  const string InputTrees = DIR["input"] + "InputTrees.txt";
  map<string, vector<string> > InputFileCollection;
  if(!getInputFileNames(InputTrees, InputFileCollection)){ return; }
  
  TObjArray* aDSTAG = new TObjArray(); // Array to store the different tags in the list of trees
  aDSTAG->SetOwner(true);
  map<string, RooWorkspace> Workspace;

  for(map<string, vector<string> >::iterator FileCollection=InputFileCollection.begin(); FileCollection!=InputFileCollection.end(); ++FileCollection) {
    // Get the file tag which has the following format: DSTAG_COLL , i.e. DATA_PP 
    string FILETAG = FileCollection->first;  
    string DSTAG   = FILETAG;
    if (FILETAG.size()) {
      DSTAG.erase(DSTAG.find("_"));
    } else {
      cout << "[ERROR] FILETAG is empty!" << endl;
    }
    // Extract the filenames
    vector<string> InputFileNames = FileCollection->second; 
    string         OutputFileName;
    // If we have data, check if the user wants to fit data
    if ( (FILETAG.find("DATA")!=std::string::npos) && fitData==true ) {
      if ( (FILETAG.find("PP")!=std::string::npos)   && !fitPP   ) continue; // If we find PP, check if the user wants PP
      if ( (FILETAG.find("PbPb")!=std::string::npos) && !fitPbPb ) continue; // If we find PbPb, check if the user wants PbPb
      OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + ".root";
      if(!tree2DataSet(Workspace[DSTAG], InputFileNames, FILETAG, OutputFileName)){ return; }
      if (!aDSTAG->FindObject(DSTAG.c_str())) aDSTAG->Add(new TObjString(DSTAG.c_str()));
    }
    // If we find MC, check if the user wants to fit MC
    if ( (FILETAG.find("MC")!=std::string::npos) && fitData==false ) {
      if ( (FILETAG.find("PP")!=std::string::npos)    && !fitPP    ) continue; // If we find PP, check if the user wants PP
      if ( (FILETAG.find("PbPb")!=std::string::npos)  && !fitPbPb  ) continue; // If we find PbPb, check if the user wants PbPb
      if ( (FILETAG.find("JPSI")!=std::string::npos)  && !incJpsi  ) continue; // If we find Jpsi MC, check if the user wants to include Jpsi
      if ( (FILETAG.find("PSI2S")!=std::string::npos) && !incPsi2S ) continue; // If we find Psi2S MC, check if the user wants to include Psi2S
      OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + ".root";
      if(!tree2DataSet(Workspace[DSTAG], InputFileNames, FILETAG, OutputFileName)){ return; }
      if (!aDSTAG->FindObject(DSTAG.c_str())) aDSTAG->Add(new TObjString(DSTAG.c_str()));
      if (wantPureSMC)
      {
        OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + "_PureS" + ".root";
        if(!tree2DataSet(Workspace[Form("%s_PureS",DSTAG.c_str())], InputFileNames, FILETAG, OutputFileName)){ return; }
      }
    } 
  }
  if (Workspace.size()==0) {
    cout << "[ERROR] No onia tree files were found matching the user's input settings!" << endl; return;
  }

  // -------------------------------------------------------------------------------
  // STEP 2: LOAD THE INITIAL PARAMETERS
  /*
    Input : List of initial parameters with format PT <tab> RAP <tab> CEN <tab> iniPar ... 
    Output: two vectors with one entry per kinematic bin filled with the cuts and initial parameters
  */
  
  string InputFile;
//.........这里部分代码省略.........

void ClusterDensityFromQA(const char* qaFile="")
{
    if ( TString(qaFile).BeginsWith("alien"))
    {
        TGrid::Connect("alien://");
    }

    TFile* f = TFile::Open(qaFile);


    std::vector<IntegrationRange> rangesLow;
    std::vector<IntegrationRange> rangesHigh;

    rangesLow.push_back(IntegrationRange{-10,10, 80, 85});
    rangesLow.push_back(IntegrationRange{-10,10, 80, 85});
    rangesLow.push_back(IntegrationRange{-10,10, 95,100});
    rangesLow.push_back(IntegrationRange{-10,10, 95,100});

    for ( int i = 0; i < 4; ++i )
    {
        IntegrationRange ref = rangesLow[i];

        TObjArray* a = static_cast<TObjArray*>(f->Get("MUON_QA/expert"));

        TH2* h = static_cast<TH2*>(a->FindObject(Form("hClusterHitMapInCh%d",i+1)));

        // derive other symetric ranges from that one
        
        std::vector<IntegrationRange> ranges;

        double ysize = ref.ymax - ref.ymin;
        double xsize = (ref.xmax - ref.xmin)/2.0;

        ranges.push_back(ref);
        ranges.push_back(IntegrationRange{ref.xmin,ref.xmax,-ref.ymax,-ref.ymin});
        ranges.push_back(IntegrationRange{ref.ymin+ysize/2.0,ref.ymax,ref.xmin,ref.xmax});

        TCanvas* c = new TCanvas(Form("Chamber%d",i+1),Form("Chamber%d",i+1));

        h->Draw("colz");

        std::cout << "CHAMBER " << i+1 << " LOW = ";

        for ( auto r : ranges )
        {
            double count = h->Integral(
                h->GetXaxis()->FindBin(r.xmin),
                h->GetXaxis()->FindBin(r.xmax),
                h->GetYaxis()->FindBin(r.ymin),
                h->GetYaxis()->FindBin(r.ymax)
                );

            std::cout << " " << count << "(" << r.Surface() << " cm^2)";
            std::vector<double> x = { r.xmin,r.xmax,r.xmax,r.xmin,r.xmin };
            std::vector<double> y = { r.ymax,r.ymax,r.ymin,r.ymin,r.ymax };

            TPolyLine* l = new TPolyLine(5,&x[0],&y[0]);
            l->SetLineColor(1);
            l->SetLineStyle(9);
            l->Draw();
        }

        std::cout << std::endl;

    }
}

void simall(Int_t nEvents = 1,
	    TObjArray& fDetList,
            Bool_t fVis=kFALSE,
            TString fMC="TGeant3",
	    TString fGenerator="mygenerator",
	    Bool_t fUserPList= kFALSE
	   )
{


  TString dir = getenv("VMCWORKDIR");
  TString simdir = dir + "/macros";

  TString sim_geomdir = dir + "/geometry";
  gSystem->Setenv("GEOMPATH",sim_geomdir.Data());

  TString sim_confdir = dir + "gconfig";
  gSystem->Setenv("CONFIG_DIR",sim_confdir.Data());

// Output files
  TString OutFile = "simout.root";
  TString ParFile = "simpar.root";


  // In general, the following parts need not be touched
  // ========================================================================

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();

  gSystem->Load("libGenVector");

  gSystem->Load("libGeoBase");
  gSystem->Load("libFairDB");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  gSystem->Load("libField");
  gSystem->Load("libGen");

  //----  Load specific libraries ---------------------------------------
  gSystem->Load("libEnsarbase");
  gSystem->Load("libEnsarGen");
  gSystem->Load("libEnsarData");
  gSystem->Load("libEnsarMyDet");
 
  // -----   Create simulation run   ----------------------------------------
  FairRunSim* run = new FairRunSim();
  run->SetName(fMC.Data());              // Transport engine
  run->SetOutputFile(OutFile.Data());          // Output file
  FairRuntimeDb* rtdb = run->GetRuntimeDb();

  //  R3B Special Physics List in G4 case
  if ( (fUserPList  == kTRUE ) &&
       (fMC.CompareTo("TGeant4")   == 0)
      ){
       run->SetUserConfig("g4Config.C");
       run->SetUserCuts("SetCuts.C");
   }


  // -----   Create media   -------------------------------------------------
  //run->SetMaterials("media_r3b.geo");       // Materials
  

  // Magnetic field map type
  // Int_t fFieldMap = 0;


  // Global Transformations
  //- Two ways for a Volume Rotation are supported
  //-- 1) Global Rotation (Euler Angles definition)
  //-- This represent the composition of : first a rotation about Z axis with
  //-- angle phi, then a rotation with theta about the rotated X axis, and
  //-- finally a rotation with psi about the new Z axis.
  Double_t phi,theta,psi;

  //-- 2) Rotation in Ref. Frame of the Volume
  //-- Rotation is Using Local Ref. Frame axis angles
  Double_t thetaX,thetaY,thetaZ;


  //- Global Translation  Lab. frame.
  Double_t tx,ty,tz;


  // -----   Create geometry --------------------------------------------

  if (fDetList.FindObject("MYDET") ) {
      //My Detector definition
//.........这里部分代码省略.........

//
// *** Configuration script for phi->KK analysis with 2010 runs ***
// 
// A configuration script for RSN package needs to define the followings:
//
// (1) decay tree of each resonance to be studied, which is needed to select
//     true pairs and to assign the right mass to all candidate daughters
// (2) cuts at all levels: single daughters, tracks, events
// (3) output objects: histograms or trees
//
Bool_t RsnConfigPhiTPC
(
   AliRsnAnalysisTask *task,
   Bool_t              isMC,
   Bool_t              isMix,
   Bool_t              useCentrality,
   AliRsnCutSet       *eventCuts
)
{
   if (!task) ::Error("RsnConfigPhiTPC", "NULL task");
   
   // we define here a suffix to differentiate names of different setups for the same resonance
   // and we define also the name of the list of tracks we want to select for the analysis
   // (if will fail if no lists with this name were added to the RsnInputHandler)
   const char *suffix     = "tpc";
   const char *listName   = "kaonTPC";
   Bool_t      useCharged =  kTRUE;
   Int_t       listID     = -1;
   
   // find the index of the corresponding list in the RsnInputHandler
   AliAnalysisManager        *mgr   = AliAnalysisManager::GetAnalysisManager();
   AliMultiInputEventHandler *multi = dynamic_cast<AliMultiInputEventHandler*>(mgr->GetInputEventHandler());
   if (multi) {
      TObjArray          *array = multi->InputEventHandlers();
      AliRsnInputHandler *rsn   = (AliRsnInputHandler*)array->FindObject("rsnInputHandler");
      if (rsn) {
         AliRsnDaughterSelector *sel = rsn->GetSelector();
         listID = sel->GetID(listName, useCharged);
      }
   }
   if (listID >= 0)
      ::Info ("RsnConfigPhiTPC.C", "Required list '%s' stays in position %d", listName, listID);
   else {
      ::Error("RsnConfigPhiTPC.C", "Required list '%s' absent in handler!", listName);
      return kFALSE;
   }
   
   // ----------------------------------------------------------------------------------------------
   // -- DEFINITIONS -------------------------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------
   
   // PAIR DEFINITIONS:
   // this contains the definition of particle species and charge for both daughters of a resonance,
   // which are used for the following purposes:
   // --> species is used to assign the mass to the daughter (e.g. for building invariant mass)
   // --> charge is used to select what tracks to use when doing the computation loops
   // When a user wants to compute a like-sign background, he must define also a pair definition
   // for each like-sign: in case of charged track decays, we need one for ++ and one for --
   // Last two arguments are necessary only in some cases (but it is not bad to well initialize them):
   // --> PDG code of resonance, which is used for selecting true pairs, when needed
   // --> nominal resonance mass, which is used for computing quantities like Y or Mt
   AliRsnPairDef *phi_kaonP_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455);
   AliRsnPairDef *phi_kaonP_kaonP = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455);
   AliRsnPairDef *phi_kaonM_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455);

   // PAIR LOOPS:
   // these are the objects which drive the computations and fill the output histograms
   // each one requires to be initialized with an AliRsnPairDef object, which provided masses,
   // last argument tells if the pair is for mixing or not (this can be also set afterwards, anyway)
   const Int_t     nPairs = 5;
   Bool_t          addPair[nPairs] = {1, 1, 1, 1, 1};
   AliRsnLoopPair *phiLoop[nPairs];
   phiLoop[0] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM"     , suffix), phi_kaonP_kaonM, kFALSE);
   phiLoop[1] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_true", suffix), phi_kaonP_kaonM, kFALSE);
   phiLoop[2] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_mix" , suffix), phi_kaonP_kaonM, kTRUE );
   phiLoop[3] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonP"     , suffix), phi_kaonP_kaonP, kFALSE);
   phiLoop[4] = new AliRsnLoopPair(Form("%s_phi_kaonM_kaonM"     , suffix), phi_kaonM_kaonM, kFALSE);

   // set additional option for true pairs
   // 1) we select only pairs coming from the same mother, which must have the right PDG code (from pairDef)
   // 2) we select only pairs decaying according to the right channel (from pairDef species+charge definitions)
   phiLoop[1]->SetOnlyTrue(kTRUE);
   phiLoop[1]->SetCheckDecay(kTRUE);
   
   // don't add true pairs if not MC
   if (!isMC) addPair[1] = 0;
   addPair[0] = !isMix;
   addPair[1] = !isMix;
   addPair[2] =  isMix;
   addPair[3] = !isMix;
   addPair[4] = !isMix;
   
   // ----------------------------------------------------------------------------------------------
   // -- COMPUTED VALUES & OUTPUTS -----------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------
   
   // All values which should be computed are defined here and passed to the computation objects,
   // since they define all that is computed bye each one, and, in case one output is a histogram
   // they define the binning and range for that value
   //
//.........这里部分代码省略.........

//_________________________________________________________________________________________
Int_t checkPullTree(TString pathTree,  TString pathNameThetaMap, TString pathNameSigmaMap,
                    TString mapSuffix, const Int_t collType /*0: pp, 1: pPb, 2: PbPb*/, const Bool_t plotPull = kTRUE,
                    const Double_t downScaleFactor = 1,
                    TString pathNameSplinesFile = "", TString prSplinesName = "",
                    TString fileNameTree = "bhess_PIDetaTree.root", TString treeName = "fTree")
{
  const Bool_t isNonPP = collType != 0;
  const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
  
  Bool_t recalculateExpecteddEdx = pathNameSplinesFile != "";
  
  TFile* f = 0x0;
	
  f = TFile::Open(Form("%s/%s", pathTree.Data(), fileNameTree.Data()));
  if (!f)  {
    std::cout << "Failed to open tree file \"" << Form("%s/%s", pathTree.Data(), fileNameTree.Data()) << "\"!" << std::endl;
    return -1;
  }
      
  // Extract the data Tree
  TTree* tree = dynamic_cast<TTree*>(f->Get(treeName.Data()));
  if (!tree) {
    std::cout << "Failed to load data tree!" << std::endl;
    return -1;
  }
  
  // Extract the splines, if desired
  TSpline3* splPr = 0x0;
  if (recalculateExpecteddEdx) {
    std::cout << "Loading splines to recalculate expected dEdx!" << std::endl << std::endl;
    
    TFile* fSpl = TFile::Open(pathNameSplinesFile.Data());
    if (!fSpl) {
      std::cout << "Failed to open spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
      return 0x0;
    }
    
    TObjArray* TPCPIDResponse = (TObjArray*)fSpl->Get("TPCPIDResponse");
    if (!TPCPIDResponse) {
      splPr = (TSpline3*)fSpl->Get(prSplinesName.Data());
      
      // If splines are in file directly, without TPCPIDResponse object, try to load them
      if (!splPr) {
        std::cout << "Failed to load object array from spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
        return 0x0;
      }
    }
    else {
      splPr = (TSpline3*)TPCPIDResponse->FindObject(prSplinesName.Data());
      
      if (!splPr) {
        std::cout << "Failed to load splines from file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
        return 0x0;
      }
    }
  }
  else
    std::cout << "Taking dEdxExpected from Tree..." << std::endl << std::endl;

  // Extract the correction maps
  TFile* fMap = TFile::Open(pathNameThetaMap.Data());
  if (!fMap)  {
    std::cout << "Failed to open thetaMap file \"" << pathNameThetaMap.Data() << "\"! Will not additionally correct data...." << std::endl;
  }

  TH2D* hMap = 0x0;
  
  if (fMap) {
    hMap = dynamic_cast<TH2D*>(fMap->Get(Form("hRefined%s", mapSuffix.Data())));
    if (!hMap) {
      std::cout << "Failed to load theta map!" << std::endl;
      return -1;
    }
  }

  TFile* fSigmaMap = TFile::Open(pathNameSigmaMap.Data());
  if (!fSigmaMap)  {
    std::cout << "Failed to open simgaMap file \"" << pathNameSigmaMap.Data() << "\"!" << std::endl;
    return -1;
  }

  TH2D* hThetaMapSigmaPar1 = dynamic_cast<TH2D*>(fSigmaMap->Get("hThetaMapSigmaPar1"));
  if (!hThetaMapSigmaPar1) {
    std::cout << "Failed to load sigma map for par 1!" << std::endl;
    return -1;
  }

  Double_t c0 = -1;
  TNamed* c0Info = dynamic_cast<TNamed*>(fSigmaMap->Get("c0"));
  if (!c0Info) {
    std::cout << "Failed to extract c0 from file with sigma map!" << std::endl;
    return -1;
  }

  TString c0String = c0Info->GetTitle();
  c0 = c0String.Atof();
  printf("Loaded parameter 0 for sigma: %f\n\n", c0);

  if (plotPull)
//.........这里部分代码省略.........

//
// Test config macro for RSN package.
// It configures:
// 1) a monitor for all tracks passing quality cuts
// 2) a monitor for all tracks passing quality + PID cuts
// 3) an unlike-sign invariant-mass + pt distribution for K+K- pairs
//
Bool_t RsnConfigPhiKaonTest
(
   AliRsnAnalysisTask *task,
   Bool_t              isMC
)
{
   if (!task) {
      ::Error("RsnConfigPhiKaonTest.C", "NULL task");
      return kFALSE;
   }
   
   const char *suffix = "test";
      
   // ----------------------------------------------------------------------------------------------
   // -- DEFINITIONS -------------------------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------
   
   // daughter definition for monitor loops
   // since it is intended to loop over all 'track' like objects (in the sense that we exclude V0s and cascades),
   // we initialize it using the constructor that requires an AliRsnDaughter::EType and a charge, but since
   // we want to loop over both charges, we set it to anything which is not '+' '-' or '0', which are tokens for
   // selecting only positive, only negative or only neutral
   AliRsnDaughterDef *tracks = new AliRsnDaughterDef(AliRsnDaughter::kTrack, 0);
   
   // definition of pair decay tree for phi resonance
   // here we *must* specify a particle species and a charge, in order to check the decay tree
   // last arguments are the PDG code and nominal mass of the resonance, which are needed when
   // one wants to select true pairs only and/or he wants to compute rapidity or Mt 
   AliRsnPairDef *pairDef = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455);
   
   // definition of loop objects:
   // (a) 1 monitor for all tracks passing quality cuts
   // (b) 1 monitor for all tracks passing quality+PID cuts
   // (c) 1 pair filled with all tracks passing same cuts as (b)
   // (d) 1 pair like (c) but for mixing
   // (e) 1 pair like (c) but with true pairs only
   // NOTE: (c) and (d) are instantiated with same settings, they will be made
   //       different after some settings done in second moment
   AliRsnLoopDaughter *loopQuality = new AliRsnLoopDaughter(Form("%s_mon_quality", suffix), 0, tracks);
   AliRsnLoopDaughter *loopPID     = new AliRsnLoopDaughter(Form("%s_mon_pid"    , suffix), 0, tracks);
   AliRsnLoopPair     *loopPhi     = new AliRsnLoopPair    (Form("%s_unlike"     , suffix), pairDef);
   AliRsnLoopPair     *loopPhiMix  = new AliRsnLoopPair    (Form("%s_unlike"     , suffix), pairDef);
   AliRsnLoopPair     *loopPhiTrue = new AliRsnLoopPair    (Form("%s_trues"      , suffix), pairDef);
   
   // set additional option for true pairs (slot [0])
   loopPhiTrue->SetOnlyTrue(kTRUE);
   loopPhiTrue->SetCheckDecay(kTRUE);
   
   // set mixing options
   loopPhi    ->SetMixed(kFALSE);
   loopPhiMix ->SetMixed(kTRUE);
   loopPhiTrue->SetMixed(kFALSE);
   
   // assign the ID of the entry lists to be used by each pair to get selected daughters
   // in our case, the AliRsnInputHandler contains only one list for selecting kaons
   Int_t idQuality, idPID;
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   AliMultiInputEventHandler *multi = dynamic_cast<AliMultiInputEventHandler*>(mgr->GetInputEventHandler());
   if (!multi) {
      myError("Needed a multi input handler!");
      return kFALSE;
   }
   TObjArray *array = multi->InputEventHandlers();
   AliRsnInputHandler *rsn = (AliRsnInputHandler*)array->FindObject("rsnInputHandler");
   if (!rsn) {
      myError("Needed an RSN event handler");
      return kFALSE;
   }
   AliRsnDaughterSelector *sel = rsn->GetSelector();
   idQuality = sel->GetID("qualityTPC", kTRUE);
   idPID = sel->GetID("kaonTPC", kTRUE);
   if (idQuality < 0 || idPID < 0) {
      myError("List problems");
      return kFALSE;
   }
   loopQuality->SetListID(idQuality);
   loopPID    ->SetListID(idPID);
   loopPhi    ->SetListID(0, idPID);
   loopPhi    ->SetListID(1, idPID);
   loopPhiMix ->SetListID(0, idPID);
   loopPhiMix ->SetListID(1, idPID);
   loopPhiTrue->SetListID(0, idPID);
   loopPhiTrue->SetListID(1, idPID);
   
   // ----------------------------------------------------------------------------------------------
   // -- EVENT CUTS --------------------------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------

   // primary vertex:
   // - 2nd argument --> |Vz| range
   // - 3rd argument --> minimum required number of contributors
   // - 4th argument --> tells if TPC stand-alone vertexes must be accepted
   // we switch on the check for pileup
//.........这里部分代码省略.........

//___________________________________________________________________
Int_t extractPtResolution(TString pathNameData, TString listName,
                          Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
                          Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/,
                          Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/)
{
  if (listName == "") {
    listName = pathNameData;
    listName.Replace(0, listName.Last('/') + 1, "");
    listName.ReplaceAll(".root", "");
  }
  
  TString pathData = pathNameData;
  pathData.Replace(pathData.Last('/'), pathData.Length(), "");
  
  TH1D* hPtResolutionFit[AliPID::kSPECIES] = {0x0, };
  TH2D* hPtResolution[AliPID::kSPECIES] = {0x0, };
  THnSparse* hPtResolutionRaw[AliPID::kSPECIES] = {0x0, };

  TFile* fileData = TFile::Open(pathNameData.Data());
  if (!fileData) {
    printf("Failed to open data file \"%s\"\n", pathNameData.Data());
    return -1;
  }
  
  TObjArray* histList = (TObjArray*)(fileData->Get(listName.Data()));
  
  if (!histList) {
    printf("Failed to load list!\n");
    return -1;
  }
  
  Double_t actualLowerCentrality = -2;
  Double_t actualUpperCentrality = -2;
  
  Double_t actualLowerJetPt = -1.;
  Double_t actualUpperJetPt = -1.;
  
  Bool_t restrictJetPtAxis = (lowerJetPt >= 0 && upperJetPt >= 0);
  const Bool_t restrictCentrality = ((lowerCentrality >= -1) && (upperCentrality >= -1));
  
  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
    const TString sparseName = Form("fPtResolution_%s", AliPID::ParticleShortName(species));
    hPtResolutionRaw[species] = (THnSparse*)histList->FindObject(sparseName.Data());
    
    if (!hPtResolutionRaw[species]) {
      printf("Failed to load THnSparse for %s: %s!\n", AliPID::ParticleShortName(species), sparseName.Data());
      return -1;
    }
    
    // Set proper errors, if not yet calculated
    if (!hPtResolutionRaw[species]->GetCalculateErrors()) {
      std::cout << "Re-calculating errors of " << hPtResolutionRaw[species]->GetName() << "..." << std::endl;
      
      hPtResolutionRaw[species]->Sumw2();
      
      Long64_t nBinsPtResolutionRaw = hPtResolutionRaw[species]->GetNbins();
      Double_t binContent = 0;
      for (Long64_t bin = 0; bin < nBinsPtResolutionRaw; bin++) {
        binContent = hPtResolutionRaw[species]->GetBinContent(bin);
        hPtResolutionRaw[species]->SetBinError(bin, TMath::Sqrt(binContent));
      }
    }
    
    // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
    const Int_t lowerCentralityBinLimit = restrictCentrality ? hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->FindBin(lowerCentrality + 0.001) 
                                                                    : -1;
    const Int_t upperCentralityBinLimit = restrictCentrality ? hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->FindBin(upperCentrality - 0.001) 
                                                                    : -2;
    
    
    if (restrictCentrality) {
      actualLowerCentrality = hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
      actualUpperCentrality = hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->GetBinLowEdge(upperCentralityBinLimit);
      
      hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
    }
    
    const Bool_t restrictCharge = (chargeMode != kAllCharged);
    
    Int_t lowerChargeBinLimit = -1;
    Int_t upperChargeBinLimit = -2;
      
    if (restrictCharge) {
      // Add subtract a very small number to avoid problems with values right on the border between to bins
      if (chargeMode == kNegCharge) {
        lowerChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(-1. + 0.001);
        upperChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(0. - 0.001);
      }
      else if (chargeMode == kPosCharge) {
        lowerChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(0. + 0.001);
        upperChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(1. - 0.001);
      }
      
      // Check if the values look reasonable
      if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1
          && upperChargeBinLimit <= hPtResolutionRaw[species]->GetAxis(kPtResCharge)->GetNbins()) {
        // OK
      }
      else {
//.........这里部分代码省略.........

void printCalibStat(Int_t run, const char * fname,  TTreeSRedirector * pcstream){

  //
  // Dump the statistical information about all histograms in the calibration files 
  //    into the statistical tree, print on the screen (log files) as well 
  //
  //
  // 1. Default dump for all histograms
  //    Information to dump:
  //    stat =Entries, Mean, MeanError,  RMS, MaxBin
  //    Branch naming convention:
  //    <detName>_<hisName><statName>
  //
  // 2. Detector statistical information  - to be implemented by expert
  //                                      - First version implemented by MI 
  //  
  // 

  TFile *fin = TFile::Open(fname);
  if (!fin) return;
  const Int_t kMaxHis=10000;
  
  TList * keyList = fin->GetListOfKeys();
  Int_t nkeys=keyList->GetEntries();
  Double_t *hisEntries = new Double_t[kMaxHis];
  Double_t *hisMean = new Double_t[kMaxHis];
  Double_t *hisMeanError = new Double_t[kMaxHis];
  Double_t *hisRMS = new Double_t[kMaxHis];
  Double_t *hisMaxBin = new Double_t[kMaxHis];
  Int_t counter=0;
  
  if (pcstream) (*pcstream)<<"calibStatAll"<<"run="<<run;
  for (Int_t ikey=0; ikey<nkeys; ikey++){
    TObject * object = fin->Get(keyList->At(ikey)->GetName());
    if (!object) continue;
    if (object->InheritsFrom("TCollection")==0) continue;
    TSeqCollection *collection  = (TSeqCollection*)object; 
    Int_t nentries= collection->GetEntries();
    for (Int_t ihis=0; ihis<nentries; ihis++){
      TObject * ohis = collection->At(ihis);
      if (!ohis) continue;
      if (ohis->InheritsFrom("TH1")==0) continue;
      TH1* phis = (TH1*)ohis;
      hisEntries[counter]=phis->GetEntries();	
      Int_t idim=1;
      if (ohis->InheritsFrom("TH2")) idim=2;
      if (ohis->InheritsFrom("TH3")) idim=3;        
      hisMean[counter]=phis->GetMean(idim);	
      hisMeanError[counter]=phis->GetMeanError(idim);	
      hisRMS[counter]=phis->GetRMS(idim);	
      hisMaxBin[counter]=phis->GetXaxis()->GetBinCenter(phis->GetMaximumBin());	
      if (pcstream) (*pcstream)<<"calibStatAll"<<
		      Form("%s_%sEntries=",keyList->At(ikey)->GetName(), phis->GetName())<<hisEntries[counter]<<	
		      Form("%s_%sMean=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMean[counter]<<	
		      Form("%s_%sMeanError=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMeanError[counter]<<	
		      Form("%s_%sRMS=",keyList->At(ikey)->GetName(), phis->GetName())<<hisRMS[counter]<<	
		      Form("%s_%sMaxBin=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMaxBin[counter];	
      //printf("Histo:\t%s_%s\t%f\t%d\n",keyList->At(ikey)->GetName(), phis->GetName(), hisEntries[counter],idim);
      counter++;
    }
    delete object;
  }    
  
  //
  // Expert dump example (MI first iteration):
  //
  // 0.)  TOF dump
  //

  Int_t tofEvents=0;
  Int_t tofTracks=0;
  TList * TOFCalib = (TList*)fin->Get("TOFHistos");      
  if (TOFCalib) {
    TH1 *histoEvents = (TH1*)TOFCalib->FindObject("hHistoVertexTimestamp");
    TH1 *histoTracks = (TH1*)TOFCalib->FindObject("hHistoDeltatTimestamp");
    if (histoEvents && histoTracks){
      tofEvents = TMath::Nint(histoEvents->GetEntries());
      tofTracks = TMath::Nint(histoTracks->GetEntries());
    }
    delete TOFCalib;
  }
  printf("Monalisa TOFevents\t%d\n",tofEvents);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFevents="<<tofEvents;
  printf("Monalisa TOFtracks\t%d\n",tofTracks);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFtracks="<<tofTracks;

  //
  // 1.)  TPC  dump - usefull events/tracks  for the calibration
  //
  Int_t tpcEvents=0;
  Int_t tpcTracks=0;
  TObject* obj = dynamic_cast<TObject*>(fin->Get("TPCCalib"));
  TObjArray* array = dynamic_cast<TObjArray*>(obj);
  TDirectory* dir = dynamic_cast<TDirectory*>(obj);
  AliTPCcalibTime  * calibTime = NULL;
  if (dir) {
    calibTime = dynamic_cast<AliTPCcalibTime*>(dir->Get("calibTime"));
  }
  else if (array){
    calibTime = (AliTPCcalibTime *)array->FindObject("calibTime");
//.........这里部分代码省略.........

void sofiaall(TString InFiles = "runlist.dat",
	    TObjArray& fDetList )
{

  // Output files
  TString outFile = "./r3bunpack.root";



  // In general, the following parts need not be touched
  // ========================================================================

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------


  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();
  gSystem->Load("libFairTools");
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libR3BMbs");
  gSystem->Load("libMbsAPI");
  gSystem->Load("libSOFIAMCStack");
  gSystem->Load("libR3BRootEvent");
  gSystem->Load("libR3BLANDEvent");
  gSystem->Load("libR3BUnpack");
  //gSystem->Load("libR3BMusicEvent");
  //gSystem->Load("libR3BProEvent");
  //gSystem->Load("libR3BCrateEvent");

 
  // -----   Create analysis run   ----------------------------------------
  FairRunAna *fRun= new FairRunAna();
  fRun->SetOutputFile(outFile);

  // ----  Load MBS
  MBSUnpack *MBSunpack= new MBSUnpack("MBS unpack", InFiles);
  fRun->AddTask(MBSunpack);

  //SOFIA Crate 1
  if (fDetList.FindObject("Crate1") ) {
    CrateUnpack *CRATEunpack= new CrateUnpack("unpack");
    fRun->AddTask(CRATEunpack);
  }

  //SOFIA Crate 2
  if (fDetList.FindObject("Crate2") ) {
    MUSICUnpack *MUSICunpack= new MUSICUnpack("MUSIC unpack");
    fRun->AddTask(MUSICunpack);
  }

  // SOFIA Tof for LCP
  if (fDetList.FindObject("LCP_TOF") ) {
    ProUnpack *PROunpack= new ProUnpack("PROTON unpack");
    fRun->AddTask(PROunpack);
  }

  // SOFIA Land Detector
  if (fDetList.FindObject("LAND") ) {
    LANDUnpack *LANDunpack= new LANDUnpack("LAND unpack");
    fRun->AddTask(LANDunpack);
  }

  // -----   Initialize analysis run   --------------------------------------
  fRun->Init();
  fRun->RunOnLmdFiles(); 
   
  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Real time " << rtime/60.0 << " min, CPU time " << ctime/60.0 
       << "min" << endl << endl;
  // ------------------------------------------------------------------------

  cout << " Test passed" << endl;
  cout << " All ok " << endl;

}

void TriggerInputsForMuonEventCuts ( TString runListFilename, TString selectedInputs="", TString defaultStorage = "raw://" )
{
  AliCDBManager::Instance()->SetDefaultStorage(defaultStorage.Data());
  TObjArray inputsList;
  inputsList.SetOwner();
  
  TObjArray* selectedInputsList = selectedInputs.Tokenize(",");

  // Read input run list
  ifstream inFile(runListFilename.Data());
  TString srun = "";
  if ( inFile.is_open() ) {
    while ( ! inFile.eof() ) {
      srun.ReadLine(inFile,kFALSE);
      if ( ! srun.IsDigit() ) continue;
      
      // For each run, read trigger inputs from OCDB
      Int_t runNumber = srun.Atoi();
      AliCDBManager::Instance()->SetRun(runNumber);
      
      // Get trigger class configuration
      AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/CTP/Config");
      if ( ! entry ) continue;
      
      THashList* runInputs = new THashList();
      runInputs->SetOwner();
      runInputs->SetUniqueID((UInt_t)runNumber);
      AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject();
      const TObjArray& trigInputsArray = trigConf->GetInputs();
      AliTriggerInput* trigInput = 0x0;
      TIter next(&trigInputsArray);
      while ( ( trigInput = static_cast<AliTriggerInput*>(next()) ) ) {
        if ( selectedInputsList->GetEntriesFast() > 0 && ! selectedInputsList->FindObject(trigInput->GetName()) ) continue;
        Int_t inputId = (Int_t)TMath::Log2(trigInput->GetMask());
        TObjString* currInput = new TObjString(trigInput->GetName());
        currInput->SetUniqueID(inputId);
        runInputs->Add(currInput);
      }
      inputsList.Add(runInputs);
    }
    inFile.close();
  }
  delete selectedInputsList;
  
  // Loop on the trigger inputs
  // and group runs with an equal list of inputs
  Int_t nentries = inputsList.GetEntries();
  TArrayI checkMask(nentries);
  checkMask.Reset(1);
  for ( Int_t irun=0; irun<nentries; irun++ ) {
    if ( checkMask[irun] == 0 ) continue;
    THashList* currList = static_cast<THashList*>(inputsList.At(irun));
    TString runRange = Form("Run range: %u", currList->GetUniqueID());
    for ( Int_t jrun=irun+1; jrun<nentries; jrun++ ) {
      if ( checkMask[jrun] == 0 ) continue;
      THashList* checkList = static_cast<THashList*>(inputsList.At(jrun));
      Bool_t isDifferent = kFALSE;
      for ( Int_t itrig=0; itrig<currList->GetEntries(); itrig++ ) {
        TObjString* currInput = static_cast<TObjString*>(currList->At(itrig));
        TObject* checkInput = checkList->FindObject(currInput->GetName());
        if ( ! checkInput || checkInput->GetUniqueID() != currInput->GetUniqueID() ) {
          isDifferent = kTRUE;
          break;
        }
      } // loop on trigger inputs
      if ( isDifferent ) continue;
      checkMask[jrun] = 0;
      runRange += Form(",%u", checkList->GetUniqueID());
    } // loop on runs
    
    TString outString = "\nSetTrigInputsMap(\"";
    for ( Int_t itrig=0; itrig<currList->GetEntries(); itrig++ ) {
      TObjString* currInput = static_cast<TObjString*>(currList->At(itrig));
      outString += Form("%s:%u,",currInput->GetString().Data(), currInput->GetUniqueID());
    }
    outString.Append("\");\n");
    outString.ReplaceAll(",\"","\"");
    outString += runRange;
    printf("%s\n", outString.Data());
  } // loop on runs
}

void extractJetTrack(){

   //TString filename = "dj_RECOPAT_18_1_rhi";
   TString filename = "dj_HCPR-GoodTrk1123_All0";
   //TString filename = "dj_HydjetQ_DJQ80_F10GSR_GoodTrk1123";
   //TString filename = "dj_HydjetQ_DJUQ80_F10GSR_GoodTrk1123";

  TFile *f = new TFile(Form("/home/sungho/sctch101/data/jettrack/%s.root",filename.Data()));

  TTree *djtree = (TTree*) f->Get("djcalo/djTree");

  bool debug = false;
  maxSampling = 10; // max number of event-by-event histogram

  int count=0;

  // variables 
  float minpt = 0.9;
  float centMin = 0, centMax = 10;
  float njet_min = 100; 
  float njeteta_max = 2.0;
  float ajet_min = 50;
  float ajeteta_max = 2.0;


  // event number of interest
  targetEvtNum.push_back(1490824);
  targetEvtNum.push_back(2084186);
  targetEvtNum.push_back(2983992);


  // prepare for output files, histograms, etc;
  prepareHist();

  // event-by-event
  Nevt = djtree->GetEntries();

  float cent = 0;
  djtree->SetBranchAddress("cent",&cent);
  
  int runNum = 0, lumNum = 0, evtNum = 0;
  djtree->SetBranchAddress("run",&runNum);
  djtree->SetBranchAddress("lumi",&lumNum);
  djtree->SetBranchAddress("evt",&evtNum);


  // tracks and jet
  int nTrk = 0; djtree->SetBranchAddress("evtnp",&nTrk);

  // each jet 
  float njeteta=0, njetphi=0, njet=0;  // near side jet
  float ajeteta=0, ajetphi=0, ajet=0;  // away side jet 

  djtree->SetBranchAddress("nljeta",&njeteta);
  djtree->SetBranchAddress("nljphi",&njetphi);
  djtree->SetBranchAddress("nljet",&njet);

  djtree->SetBranchAddress("aljeta",&ajeteta);
  djtree->SetBranchAddress("aljphi",&ajetphi);
  djtree->SetBranchAddress("aljet",&ajet);

  
  // each tracks
  djtree->SetBranchAddress("ppt",&trkpt);
  djtree->SetBranchAddress("peta",&trketa);
  djtree->SetBranchAddress("pphi",&trkphi);


  for(Long_t i=0;i<djtree->GetEntries();i++){

    djtree->GetEntry(i);
    if((i%100)==0) cout<<"counting every 100 events = "<<i<<endl;

    if(debug) cout<<"Evt: "<<evtNum<<" Near side jet Et = "<<njet<<" number of tracks = "<<nTrk<<endl;

    if(cent<centMin || cent>centMax) continue; // centrality 
    if(njet<njet_min || njet>500) continue; // near side jet et cut
    if(fabs(njeteta)>njeteta_max) continue; // near side jet eta cut
    if(ajet<ajet_min) continue;
    if(fabs(ajeteta)>ajeteta_max) continue;
    
    float dphi = nljphi - aljphi;
    if(fabs(dphi)>=(TMath::Pi())) dphi = 2.*TMath::Pi() - fabs(dphi); 
    if(dphi>(TMath::Pi()*(5/6)) continue;  // dphi cut for back-to-back jets

    for(Long_t j=0;j<2;j++){
      
      float jet = (j==0) ? njet : ajet;
      float jeta = (j==0) ? njeteta : ajeteta;
      float jphi = (j==0) ? njetphi : ajetphi;

      // randomize
      //float jeta = (j==0) ? rdn.Uniform(-2,2) : rdn.Uniform(-2,2);
      //float jphi = 0;
      //if(j==0) jphi = rdn.Uniform(0,TMath::Pi());
      //else jphi = -1.*jphi; //back to back
      
      if(hdEtadPhiJetArray.FindObject(Form("hdEtadPhiJet_%d",count)) && count<maxSampling)
	((TH2F*) hdEtadPhiJetArray.FindObject(Form("hdEtadPhiJet_%d",count)))->Fill(jeta,jphi);

//.........这里部分代码省略.........

//get the RooDataSet, multiply each roorealvar by the event_weight
void get_roodset_from_ttree(TDirectoryFile *f, TString treename, RooDataSet* &roodset){
  cout << "Creating roodset from file: " << f->GetName() << " with tree: " << treename.Data() << endl;

  TTree *t = NULL;
  assert(roodset==NULL);
  f->GetObject(treename.Data(),t);
  if (!t) {cout << "Impossible to find TTree " << treename.Data() << endl; return;}
  TObjArray *objs = t->GetListOfBranches();
  //disables all branches
  t->SetBranchStatus("*",0);


  float v_rooisopv1;
  float v_rooisopv2;
  float v_rooisowv1;
  float v_rooisowv2;
  float v_roovar1;
  float v_roovar2;
  float v_roopt1;
  float v_roosieie1;
  float v_rooeta1;
  float v_roopt2;
  float v_roosieie2;
  float v_rooeta2;
  float v_roodiphopt;
  float v_roodiphomass;
  float v_roorho;
  float v_roosigma;
  float v_roonvtx;
  float v_rooweight;

  TBranch *b_roovar1;
  TBranch *b_roovar2;
  TBranch *b_rooisopv1;
  TBranch *b_rooisopv2;
  TBranch *b_rooisowv1;
  TBranch *b_rooisowv2;
  TBranch *b_roopt1;
  TBranch *b_roosieie1;
  TBranch *b_rooeta1;
  TBranch *b_roopt2;
  TBranch *b_roosieie2;
  TBranch *b_rooeta2;
  TBranch *b_roodiphopt;
  TBranch *b_roodiphomass;
  TBranch *b_roorho;
  TBranch *b_roosigma;
  TBranch *b_roonvtx;
  TBranch *b_rooweight;

  const int nvars = 18;
  float* ptrs[nvars]={&v_roovar1,&v_roovar2,&v_rooisopv1,&v_rooisopv2,&v_rooisowv1,&v_rooisowv2,&v_roopt1,&v_roosieie1,&v_rooeta1,&v_roopt2,&v_roosieie2,&v_rooeta2,&v_roodiphopt,&v_roodiphomass,&v_roorho,&v_roosigma,&v_roonvtx,&v_rooweight};
  TBranch** branches[nvars]={&b_roovar1,&b_roovar2,&b_rooisopv1,&b_rooisopv2,&b_rooisowv1,&b_rooisowv2,&b_roopt1,&b_roosieie1,&b_rooeta1,&b_roopt2,&b_roosieie2,&b_rooeta2,&b_roodiphopt,&b_roodiphomass,&b_roorho,&b_roosigma,&b_roonvtx,&b_rooweight};
  RooRealVar* rooptrs[nvars]={roovar1,roovar2,rooisopv1,rooisopv2,rooisowv1,rooisowv2,roopt1,roosieie1,rooeta1,roopt2,roosieie2,rooeta2,roodiphopt,roodiphomass,roorho,roosigma,roonvtx,rooweight};
  bool status[nvars];
  RooArgSet args;
  for (int i=0; i<nvars; i++){ 
    status[i]=0;
    TString name = rooptrs[i]->GetName();
    TObject *obj = objs->FindObject(name.Data());
    if (!obj) continue;
    t->SetBranchStatus(name.Data(),1);
    status[i]=1;
    t->SetBranchAddress(name.Data(),ptrs[i],branches[i]);
 
   args.add(*(rooptrs[i]));
  }

  TString newname = Form("roo_%s",t->GetName());
  roodset = new RooDataSet(newname.Data(),newname.Data(),args,WeightVar(*rooweight) );
  for (int j=0; j<t->GetEntries(); j++){
    t->GetEntry(j);
    for (int i=0; i<nvars; i++){
      if (!status[i]) continue;
      rooptrs[i]->setVal(*(ptrs[i]));
    }
    roodset->add(args,v_rooweight);
  }

  cout << "Imported roodset " << newname.Data() << " from TTree " << t->GetName() << endl;
  roodset->Print();

}