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

cmsListDir( const TObjString * firstdirname, const TDirectory * firstDir )
 {
  TObjArray * dirs = new TObjArray ;
  dirs->AddLast(new TPair(firstdirname,firstDir)) ;
  TList * keys ;
  TKey * key ;
  TH1 * histo ;
  TIter nextDir(dirs) ;
  TPair * pair ;
  const TObjString * dirname ;
  const TDirectory * dir ;
  while (pair = (TPair *)nextDir())
   {
    dirname = (TObjString *)pair->Key() ;
    dir = (TDirectory *)pair->Value() ;
    keys = dir->GetListOfKeys() ;
    TIter nextKey(keys) ;
    while (key = (TKey *)nextKey())
     {
      obj = key->ReadObj() ;
      if (obj->IsA()->InheritsFrom("TDirectory"))
       {
        dirs->AddLast(new TPair(new TObjString(dirname->String()+"/"+obj->GetName()),obj)) ;
       }
      else if (obj->IsA()->InheritsFrom("TH1"))
       {
        histo = (TH1 *)obj ;
        std::cout
          <<"Histo "<<dirname->String()<<"/"<<histo->GetName()
          <<" has "<<histo->GetEffectiveEntries()<<" entries"
          <<" of mean value "<<histo->GetMean()
          <<std::endl ;
       }
      else
       { std::cout<<"What is "<<obj->GetName()<<" ?"<<std::endl ; }
     }
   }
 }

//____________________________________________________
TObjArray* CreateRecoParamObject()
{
  TObjArray *recos = new TObjArray(4);

  AliTRDrecoParam *rec = 0x0;
  recos->AddLast(rec = AliTRDrecoParam::GetLowFluxParam());
  rec->SetEventSpecie(AliRecoParam::kLowMult);
  rec->SetNameTitle("LOW", "TRD Low Flux Reco Param");
  rec->SetStreamLevel(AliTRDrecoParam::kTracker, 1);
  rec->SetAsDefault();
  // further settings for low flux reco param
  // reco->SetThisAndThat()

  recos->AddLast(rec = AliTRDrecoParam::GetHighFluxParam());
  rec->SetEventSpecie(AliRecoParam::kHighMult);
  rec->SetNameTitle("HIGH", "TRD High Flux Reco Param");
  rec->SetStreamLevel(AliTRDrecoParam::kTracker, 1);

  recos->AddLast(rec = AliTRDrecoParam::GetCosmicTestParam());
  rec->SetEventSpecie(AliRecoParam::kCosmic);
  rec->SetNameTitle("COSMIC", "TRD Cosmic Reco Param");
  rec->SetStreamLevel(AliTRDrecoParam::kTracker, 1);
  rec->SetRawStreamVersion("FAST");

  recos->AddLast(rec = AliTRDrecoParam::GetCosmicTestParam());
  rec->SetEventSpecie(AliRecoParam::kCalib);
  rec->SetNameTitle("CALIBRATION", "TRD Calibration Reco Param");
  rec->SetStreamLevel(AliTRDrecoParam::kTracker, 1);
  rec->SetRawStreamVersion("FAST");

//  recos->AddLast(rec = AliTRDrecoParam::GetLowFluxParam());
//  rec->SetNameTitle("HLT", "TRD HLT Reco Param");
//  rec->SetChi2Y(.1);
//  rec->SetChi2Z(5.);

  return recos;
}

void ilceve_init_import_macros()
{
  // Put macros in the list of browsables, add a macro browser to
  // top-level GUI.

  TString macdir("$(ILC_ROOT)/EVE/ilc-macros");
  gSystem->ExpandPathName(macdir);

  TFolder* f = gEve->GetMacroFolder();
  void* dirhandle = gSystem->OpenDirectory(macdir.Data());
  if (dirhandle != 0)
  {
    char* filename;
    TPMERegexp re("\\.C$");
    TObjArray names;
    while ((filename = gSystem->GetDirEntry(dirhandle)) != 0)
    {
      if (re.Match(filename))
	names.AddLast(new TObjString(filename));
    }
    names.Sort();

    for (Int_t ii=0; ii<names.GetEntries(); ++ii)
    {
      TObjString * si = (TObjString*) names.At(ii);
       f->Add(new TEveMacro(Form("%s/%s", macdir.Data(), (si->GetString()).Data())));
    }
  }
  gSystem->FreeDirectory(dirhandle);

  gROOT->GetListOfBrowsables()->Add(new TSystemDirectory(macdir.Data(), macdir.Data()));

  {
    TEveBrowser   *br = gEve->GetBrowser();
    TGFileBrowser *fb = 0;
    fb = br->GetFileBrowser();
    fb->GotoDir(macdir);
    {
      br->StartEmbedding(0);
      fb = br->MakeFileBrowser();
      fb->BrowseObj(f);
      fb->Show();
      br->StopEmbedding();
      br->SetTabTitle("Macros", 0);
      br->SetTab(0, 0);
    }
  }
}

void Example_tags(TString topDir = "/star/rcf/GC/daq/tags")
{
//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Example_tags.C                                                       //
//                                                                      //
// shows how to use the STAR tags files                                 //
// Input: top level directory                                           //
//                                                                      //
// what it does:                                                        //
// 1. creates TChain from all tags files down from the topDir           //
// 2. loops over all events in the chain                                //
//                                                                      //
// owner: Alexandre V. Vaniachine <[email protected]>                //
//////////////////////////////////////////////////////////////////////////

  gSystem->Load("libTable");
  gSystem->Load("St_base");
  // start benchmarks
  gBenchmark = new TBenchmark();
  gBenchmark->Start("total");
   
  // set loop optimization level
  gROOT->ProcessLine(".O4");
  // gather all files from the same top directory into one chain
  // topDir must end with "/"
  topDir +='/';
  St_FileSet dirs(topDir);
  St_DataSetIter next(&dirs,0);
  St_DataSet *set = 0; 
  TChain chain("Tag");
  while ( (set = next()) ) {           
    if (strcmp(set->GetTitle(),"file") || 
	!(strstr(set->GetName(),".tags.root"))) continue;
    chain.Add(gSystem->ConcatFileName(topDir,set->Path()));
  }
  UInt_t nEvents = chain->GetEntries();
  cout<<"chained "<<nEvents<<" events "<<endl;

  TObjArray *files = chain.GetListOfFiles();
  UInt_t nFiles = files->GetEntriesFast();
  cout << "chained " << nFiles << " files from " << topDir << endl;

  TObjArray *leaves = chain.GetListOfLeaves();
  Int_t nleaves = leaves->GetEntriesFast();

  TString tableName = " ";
  TObjArray *tagTable = new TObjArray;

  Int_t tableCount = 0;
  Int_t *tableIndex = new Int_t[nleaves];
  Int_t tagCount = 0;

  // decode tag table names
  for (Int_t l=0;l<nleaves;l++) {
    TLeaf *leaf = (TLeaf*)leaves->UncheckedAt(l);
    tagCount+=leaf->GetNdata();
    TBranch *branch = leaf->GetBranch();
    // new tag table name
    if ( strstr(branch->GetName(), tableName.Data()) == 0 ) {
      tableName = branch->GetName();
      // the tableName is encoded in the branch Name before the "."
      tableName.Resize(tableName->Last('.'));
      tagTable->AddLast(new TObjString(tableName.Data()));
      tableCount++;
    }
    tableIndex[l]=tableCount-1;
  }
  cout << " tot num tables, tags = " << tableCount << "   " 
       << tagCount << endl << endl;

  //EXAMPLE 1: how to print out names of all tags and values for first event
  for (l=0;l<nleaves;l++) {
    leaf = (TLeaf*)leaves->UncheckedAt(l);
    branch = leaf->GetBranch();
    branch->GetEntry();
    // tag comment is in the title
    TString Title = leaf->GetTitle();
    Int_t dim = leaf->GetNdata();
    if (dim==1) {
      Title.ReplaceAll('['," '"); 
      Title.ReplaceAll(']',"'"); 
    }
    cout << "\n Table: ";
    cout.width(10);
    cout << ((TObjString*)tagTable->UncheckedAt(tableIndex[l]))->GetString()
	 <<" -- has tag: " << Title << endl;
    for (Int_t i=0;i<dim;i++) {
      cout <<"                               "<< leaf->GetName();
      if (dim>1) cout << '['<<i<<']';
      cout << " = " << leaf->GetValue(i) << endl; 
    }
  }

  // EXAMPLE 2: how to make a plot
  c1 = new TCanvas("c1","Beam-Gas Rejection",600,1000);
  gStyle->SetMarkerStyle(8);
  chain->Draw("n_trk_tpc[0]:n_trk_tpc[1]");

  // EXAMPLE 3: how to make a selection (write selected event numbers on the plot)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	container -> SetStepTitle(8, "RecoCuts");
	container -> SetStepTitle(9, "RecoPID");


	//CREATE THE  CUTS -----------------------------------------------

	// Gen-Level kinematic cuts
	AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");

	//Particle-Level cuts:
	AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
	Bool_t useAbsolute = kTRUE;
	if (isSign != 2){
		useAbsolute = kFALSE;
	}
	mcGenCuts->SetRequirePdgCode(pdgCode, useAbsolute);  // kTRUE set in order to include antiparticle
	mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)

	// Acceptance cuts:
	AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts");
	AliCFTrackKineCuts *kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts");
	kineAccCuts->SetPtRange(ptmin,ptmax);
	kineAccCuts->SetEtaRange(etamin,etamax);

	// Rec-Level kinematic cuts
	AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");

	AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");

	AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");

	printf("CREATE MC KINE CUTS\n");
	TObjArray* mcList = new TObjArray(0) ;
	mcList->AddLast(mcKineCuts);
	mcList->AddLast(mcGenCuts);

	printf("CREATE ACCEPTANCE CUTS\n");
	TObjArray* accList = new TObjArray(0) ;
	accList->AddLast(kineAccCuts);

	printf("CREATE RECONSTRUCTION CUTS\n");
	TObjArray* recList = new TObjArray(0) ;   // not used!!
	recList->AddLast(recKineCuts);
	recList->AddLast(recQualityCuts);
	recList->AddLast(recIsPrimaryCuts);

	TObjArray* emptyList = new TObjArray(0);

	//CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
	printf("CREATE INTERFACE AND CUTS\n");
	AliCFManager* man = new AliCFManager() ;
	man->SetParticleContainer(container);
	man->SetParticleCutsList(0 , mcList); // MC, Limited Acceptance
	man->SetParticleCutsList(1 , mcList); // MC
	man->SetParticleCutsList(2 , accList); // Acceptance
	man->SetParticleCutsList(3 , emptyList); // Vertex
	man->SetParticleCutsList(4 , emptyList); // Refit
	man->SetParticleCutsList(5 , emptyList); // AOD
	man->SetParticleCutsList(6 , emptyList); // AOD in Acceptance
	man->SetParticleCutsList(7 , emptyList); // AOD with required n. of ITS clusters
	man->SetParticleCutsList(8 , emptyList); // AOD Reco (PPR cuts implemented in Task)
	man->SetParticleCutsList(9 , emptyList); // AOD Reco PID

	// Get the pointer to the existing analysis manager via the static access method.
	//==============================================================================
	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();

//.........这里部分代码省略.........
            double EfficiencyBinValue = STABinValue / ParticleBinValue * 100.;
            double EfficiencyBinError = sqrt(EfficiencyBinValue * (100. - EfficiencyBinValue) / ParticleBinValue);
            cout << ParticleBinValue << ", " << STABinValue << ", " << EfficiencyBinValue << endl;
            Efficiency2simPtHist->SetBinContent(PtIndex, EfficiencyBinValue);
            Efficiency2simPtHist->SetBinError(PtIndex, EfficiencyBinError);

            TH1D* ChargeCheckHist = ChargeCheck2simPtHist->ProjectionY("ChargeCheck", PtIndex, PtIndex, "o");
            double ReverseChargeBinValue = ChargeCheckHist->GetBinContent(2);
            double CoverseChargeBinValue = ChargeCheckHist->GetBinContent(4);
            double TotalChargeBinValue = ReverseChargeBinValue + CoverseChargeBinValue;
            if(TotalChargeBinValue == 0.);
                TotalChargeBinValue += 1.;
            double ReverseChargeRato = ReverseChargeBinValue / TotalChargeBinValue;
            InverseChargeRato2simPtHist->SetBinContent(PtIndex, ReverseChargeRato);

            TH1D* DeltaPtHist = DeltaPtmaxPurity2simPtHist->ProjectionY("DeltaPt", PtIndex, PtIndex, "o");
            double DeltaPtMean = DeltaPtHist->GetMean();
            double DeltaPtRMS = DeltaPtHist->GetRMS();
            DeltaPt2simPtHist->SetBinContent(PtIndex, DeltaPtMean);
            DeltaPt2simPtHist->SetBinError(PtIndex, DeltaPtRMS);

            TH1D* DeltaPhiHist = DeltaPhimaxPurity2simPtHist->ProjectionY("DeltaPhi", PtIndex, PtIndex, "o");
            double DeltaPhiMean = DeltaPhiHist->GetMean();
            double DeltaPhiRMS = DeltaPhiHist->GetRMS();
            DeltaPhi2simPtHist->SetBinContent(PtIndex, DeltaPhiMean);
            DeltaPhi2simPtHist->SetBinError(PtIndex, DeltaPhiRMS);

            TH1D* DeltaEtaHist = DeltaEtamaxPurity2simPtHist->ProjectionY("DeltaEta", PtIndex, PtIndex, "o");
            double DeltaEtaMean = DeltaEtaHist->GetMean();
            double DeltaEtaRMS = DeltaEtaHist->GetRMS();
            DeltaEta2simPtHist->SetBinContent(PtIndex, DeltaEtaMean);
            DeltaEta2simPtHist->SetBinError(PtIndex, DeltaEtaRMS);
        }
        myEfficiencyHist->AddLast(Efficiency2simPtHist);
        myParticleHist->AddLast(Particle2simPtHist);
        mySTAHist->AddLast(STA2simPtHist);
        myChargeCheckHist->AddLast(InverseChargeRato2simPtHist);
        myDeltaPtHist->AddLast(DeltaPt2simPtHist);
    }
    double minX = 0;
    double minY = 0;
    double maxX = 110;
    double maxY = 40;

    TCanvas* myCanvas = new TCanvas("Canvas", "Canvas", 800, 600);
    myCanvas->cd();
    TPad* myPad = new TPad("Pad", "Pad", 0, 0, 1, 1);
    myPad->Draw();
    myPad->cd();

    ((TH1D*)(myParticleHist->At(0)))->SetStats(0);
    ((TH1D*)(myParticleHist->At(0)))->GetXaxis()->SetTitle("simPt/Gev");
    ((TH1D*)(myParticleHist->At(0)))->GetXaxis()->CenterTitle(1);
    ((TH1D*)(myParticleHist->At(0)))->Draw();
    for(int Index = 0; Index < FileNumber; Index++) {
	((TH1D*)(mySTAHist->At(Index)))->SetStats(0);
        ((TH1D*)(mySTAHist->At(Index)))->SetLineColor(kRed+Index);
        ((TH1D*)(mySTAHist->At(Index)))->Draw("same");
    }
    TLegend *STALeg = new TLegend(0.6,0.1,0.9,0.3);
    STALeg->SetBorderSize(1);
    TString LegKey = "ParticleTrack";
    STALeg->AddEntry(myParticleHist->At(0), LegKey, "lpf");
    for(int Index = 0; Index < FileNumber; Index++) {
        LegKey = TypeName[Index];
        STALeg->AddEntry(mySTAHist->At(Index), LegKey, "lpf");

void TBDataParser::OnStartElement(const char *element, const TList *attributes)
{
	TXMLAttr *attr;
	TIter next(attributes);

	char *name = element;
	_currentElement = new TString(element);
	
	char *method = NULL;

	while ((attr = (TXMLAttr*) next())) {
		char *attrName = attr->GetName();
		char *attrValue = attr->GetValue();

		if(!strcmp(attrName, "name")) {
			name = attrValue; 
		}
		
		if(!strcmp(attrName, "method")) {
			method = attrValue;	
			_currentMethod = new TString(method);
		}
		
	}

	TFolder *currentFolder = _foldersStack->Last();

	if(!strcmp(element, "vector")) {
		TString *nameString = new TString(name);
		vector<float> *values = new vector<float>;
		_values = values;

		nameString->ReplaceAll(" ","_");

		TObjString *currentVector = new TObjString(nameString->Data());

		TObjArray *vectorsStack = _vectorsStack;
		vectorsStack->AddLast(currentVector);

		TString *joinedName = new TString(((TObjString *) vectorsStack->First())->GetString().Data());

		for(Int_t i = 1; i < vectorsStack->GetEntries(); i++) {
			joinedName->Append("_");
			joinedName->Append(((TObjString *) vectorsStack->At(i))->GetString().Data());
		}
		
		TObjString *joinedNameObj = new TObjString(joinedName->Data());
		TObjArray *joinedNameStack = _joinedNameStack;
		
		TNtuple *vector = new TNtuple(joinedName->Data(),joinedName->Data(),"values");
		_vector = vector;
		currentFolder->Add(vector);

		if(joinedNameStack->Contains(joinedName->Data()))
			cout << joinedName->Data() << "Error: " << "vector already exists. Be sure that in your XML file every vector has a unique path + name combination" << endl;
			
		joinedNameStack->Add(joinedNameObj);
			
			
		return;
	}

	_foldersStack->AddLast(currentFolder->AddFolder(name, name));
}

void getHistosFromRE(const string&   mhid,
		     const string&   filepath,
		     const string&   sre,
		     vector<std::pair<string,wTH1*> >&  v_wth1)
{
  if (gl_verbose)
    cout<<"Searching for regexp "<<sre<<" in "<<filepath;

  // allow for multiple regexes in OR combination
  //
  vector<string> v_regexes;
  Tokenize(sre,v_regexes,"|");
  if (!v_regexes.size())
    v_regexes.push_back(sre);

  // Build validated TRegexp arguments in preparation for directory recursion
  //
  TObjArray *Args = new TObjArray();
  for (size_t i=0; i<v_regexes.size(); i++) {
    TRegexp re(v_regexes[i].c_str(),kTRUE);
    if (re.Status() != TRegexp::kOK) {
      cerr << "The regexp " << v_regexes[i] << " is invalid, Status() = ";
      cerr << re.Status() << endl;
      exit(-1);
    }
    else {
      Args->AddLast(new TObjString(v_regexes[i].c_str()));
    }
  }

  // Get the root file
  //
  TFile *rootfile = openRootFile(filepath);

  if (!rootfile) {
    cerr << "File failed to open, " << filepath << endl;
    Args->Delete();
    delete Args;
    return;
  }

  // Do the recursion, collect matches
  //
  TObjArray *Matches = new TObjArray();
  recurseDirs(rootfile, &regexMatchHisto, Args, Matches);
  Args->Delete();
  delete Args;

  // Returns two objects per match: 
  // 1. the (string) path that was matched and
  // 2. the object whose path matched
  //
  int nx2matches = Matches->GetEntriesFast();
  if (gl_verbose) cout << "... " << nx2matches/2 << " match(es) found.";

  // Add the matches to the global map of histos
  //
  int istart = v_wth1.size();

  for (int i=0; i<nx2matches; i+=2) {
    TString fullspec = ((TObjString *)(*Matches)[i])->GetString();
    wTH1 *wth1 = new wTH1((TH1 *)((*Matches)[i+1]));
    wth1->histo()->UseCurrentStyle();
    wth1->histo()->SetLineColor(((i/2)%9)+1);
    wth1->histo()->SetLineStyle((i/18)+1);
    wth1->histo()->SetLineWidth(2);
    wth1->SetLegendEntry(wth1->histo()->GetName());
    string hidi= mhid+"_"+int2str(istart+(i/2));
    v_wth1.push_back(std::pair<string,wTH1 *>(hidi,wth1));

    //glmap_objpath2id.insert(pair<string,string>(fullspec,hidi));
    glmap_id2histo.insert(pair<string,wTH1 *>(hidi,wth1));
    glmap_id2objpath.insert(pair<string,string>(hidi,string(fullspec.Data())));
  }

  //Matches->Delete(); // need the histos!
  delete Matches;

  if (gl_verbose) cout << endl;
}                                                     // getHistosFromRE

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

        if(debug) cout << "SimTrackId: " << SimTrackId << ", SimTrackType: " << SimTrackType << ", PassSegmentFilter: " << PassSegmentFilter << ", SeedNumber: " << SeedNumber << ", SeedPurity: " << SeedPurity << ", SeedCharge: " << SeedCharge << ", SimTrackCharge: " << SimTrackCharge << ", SimBendingPhi: " << SimBendingPhi << ", RecBendingPhi: " << RecBendingPhi << ", RecBendingLastPhi: " << RecBendingLastPhi << ", RecMomentumPtatRef: " << RecMomentumPtatRef << ", SimMomentumPtatRef: " << SimMomentumPtatRef << ", LastSeedNumber: " << LastSeedNumber << endl;

        //if(SimTrackMomentum > 20.0)
            //continue;

        if(SeedNumber != 0) {
            if(SimTrackvalid == 1 && SeedPurity == 1. && PassSegmentFilter == true) {
                SeedPtforSimTrackvalidHist->Fill(RecMomentumPtatRef);
                SeeddeltaPtforSimTrackvalidHist->Fill((RecMomentumPtatRef-SimMomentumPtatRef)/SimMomentumPtatRef);
                SeeddeltaPhiforSimTrackvalidHist->Fill(RecDirectionPhiatRef-SimDirectionPhiatRef);
                SeeddeltaEtaforSimTrackvalidHist->Fill(RecDirectionEtaatRef-SimDirectionEtaatRef);
                RecBendingPhi2PtHist->Fill(SimMomentumPtatRef*SimTrackCharge, RecBendingPhi);
                double PtRato = SimMomentumPtatRef / RecMomentumPtatRef;
                if(debug) cout << "PtRato: " << PtRato << ", at RecBendingPhi: " << RecBendingPhi << endl;
                PtRatoofRecBendingPhiHist->Fill(RecBendingPhi, PtRato);
                RecBendingLastPhiHist->Fill((RecBendingLastPhi == 0. ? 0 : RecBendingLastPhi/fabs(RecBendingLastPhi))*SimTrackCharge);
                SeedPurityforSimTrackvalidHist->Fill(SeedPurity);
                ChargeCheckforSimTrackvalidHist->Fill(SeedCharge*SimTrackCharge);
                /*
                if(SeedPurity == 1) {
                    //SeedPtforSimTrackvalidHist->Fill(RecMomentumPtatRef);
                    //SeeddeltaPtforSimTrackvalidHist->Fill(RecMomentumPtatRef-SimMomentumPtatRef);
                    RecBendingPhi2PtHist->Fill(SimMomentumPtatRef*SimTrackCharge, RecBendingPhi);
                    double PtRato = SimMomentumPtatRef / RecMomentumPtatRef;
                    if(debug) cout << "PtRato: " << PtRato << ", at RecBendingPhi: " << RecBendingPhi << endl;
                    PtRatoofRecBendingPhiHist->Fill(RecBendingPhi, PtRato);
                    RecBendingLastPhiHist->Fill((RecBendingLastPhi == 0. ? 0 : RecBendingLastPhi/fabs(RecBendingLastPhi))*SimTrackCharge);
                }
                */
                if(SeedCharge*SimTrackCharge == -1) {
                    if(debug) cout << "R1: " << sqrt(SimBendingEntryPositionX*SimBendingEntryPositionX+SimBendingEntryPositionY*SimBendingEntryPositionY) << ", R2: " << sqrt(SimBendingLeavePositionX*SimBendingLeavePositionX+SimBendingLeavePositionY*SimBendingLeavePositionY) << endl;
                    TLine* SimSegment = new TLine(SimBendingEntryPositionX, SimBendingEntryPositionY, SimBendingLeavePositionX, SimBendingLeavePositionY);
                    SimReverseBending->AddLast(SimSegment);
                }
                LastPurityFull = true;
            }
        }
        else {
            if(LastSeedNumber != -1) {
                if(LastSimTrackvalid == 1 && LastPassSegmentFilter == true) {
                    if(debug) cout << "Filling valid track efficiency " << LastSeedNumber << endl;
                    SeedNumberforSimTrackvalidHist->Fill(LastSeedNumber);
                    SeedEfficiencyforSimTrackvalidHist->Fill(LastPurityFull==true?1:0);
                }
                else {
                    SeedNumberforSimTrackinvalidHist->Fill(LastSeedNumber);
                    SeedEfficiencyforSimTrackinvalidHist->Fill(LastPurityFull==true?1:0);
                }
                LastPurityFull = false;
            }
        }
        LastPassSegmentFilter = PassSegmentFilter;
        LastSeedNumber = SeedNumber;
        LastSimTrackvalid = SimTrackvalid;
    }
    if(LastSeedNumber != -1) {
        if(LastSimTrackvalid == 1 && LastPassSegmentFilter == true) {
            SeedNumberforSimTrackvalidHist->Fill(LastSeedNumber);
            SeedEfficiencyforSimTrackvalidHist->Fill(LastPurityFull==true?1:0);
        }
        else {
            SeedNumberforSimTrackinvalidHist->Fill(LastSeedNumber);
            SeedEfficiencyforSimTrackinvalidHist->Fill(LastPurityFull == true?1:0);
        }
    }

void limit() {
//This program demonstrates the computation of 95 % C.L. limits.
//It uses a set of randomly created histograms.
//
//Author:  [email protected] on 21.08.02

// Create a new canvas.
  TCanvas *c1 = new TCanvas("c1","Dynamic Filling Example",200,10,700,500);
  c1->SetFillColor(42);
  
// Create some histograms
  TH1D* background = new TH1D("background","The expected background",30,-4,4);
  TH1D* signal     = new TH1D("signal","the expected signal",30,-4,4);
  TH1D* data       = new TH1D("data","some fake data points",30,-4,4);
  background->SetFillColor(48);
  signal->SetFillColor(41);
  data->SetMarkerStyle(21);
  data->SetMarkerColor(kBlue);
  background->Sumw2(); // needed for stat uncertainty
  signal->Sumw2(); // needed for stat uncertainty
  
// Fill histograms randomly
  TRandom2 r;
  Float_t bg,sig,dt;
  for (Int_t i = 0; i < 25000; i++) {
     bg  = r.Gaus(0,1);
     sig = r.Gaus(1,.2);
     background->Fill(bg,0.02);
     signal->Fill(sig,0.001);
  }
  for (Int_t i = 0; i < 500; i++) {
     dt = r.Gaus(0,1);
     data->Fill(dt);
  }
  THStack *hs = new THStack("hs","Signal and background compared to data...");
  hs->Add(background);
  hs->Add(signal);
  hs->Draw("hist");
  data->Draw("PE1,Same");
  c1->Modified();
  c1->Update();
  c1->GetFrame()->SetFillColor(21);
  c1->GetFrame()->SetBorderSize(6);
  c1->GetFrame()->SetBorderMode(-1);
  c1->Modified();
  c1->Update();
  gSystem->ProcessEvents();

// Compute the limits
  cout << "Computing limits... " << endl;
  TLimitDataSource* mydatasource = new TLimitDataSource(signal,background,data);
  TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
  cout << "CLs    : "   << myconfidence->CLs()  << endl;
  cout << "CLsb   : "   << myconfidence->CLsb() << endl;
  cout << "CLb    : "   << myconfidence->CLb()  << endl;
  cout << "< CLs >  : " << myconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << myconfidence->GetExpectedCLb_b()  << endl;

// Add stat uncertainty
  cout << endl << "Computing limits with stat systematics... " << endl;
  TConfidenceLevel *mystatconfidence = TLimit::ComputeLimit(mydatasource,50000,true);
  cout << "CLs    : "   << mystatconfidence->CLs()  << endl;
  cout << "CLsb   : "   << mystatconfidence->CLsb() << endl;
  cout << "CLb    : "   << mystatconfidence->CLb()  << endl;
  cout << "< CLs >  : " << mystatconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << mystatconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << mystatconfidence->GetExpectedCLb_b()  << endl;

// Add some systematics
  cout << endl << "Computing limits with systematics... " << endl;
  TVectorD errorb(2);
  TVectorD errors(2);
  TObjArray* names = new TObjArray();
  TObjString name1("bg uncertainty");
  TObjString name2("sig uncertainty");
  names->AddLast(&name1);
  names->AddLast(&name2);
  errorb[0]=0.05; // error source 1: 5%
  errorb[1]=0;    // error source 2: 0%
  errors[0]=0;    // error source 1: 0%
  errors[1]=0.01; // error source 2: 1%
  TLimitDataSource* mynewdatasource  = new TLimitDataSource();
  mynewdatasource->AddChannel(signal,background,data,&errors,&errorb,names);
  TConfidenceLevel *mynewconfidence = TLimit::ComputeLimit(mynewdatasource,50000,true);
  cout << "CLs    : " << mynewconfidence->CLs()  << endl;
  cout << "CLsb   : " << mynewconfidence->CLsb() << endl;
  cout << "CLb    : " << mynewconfidence->CLb()  << endl;
  cout << "< CLs >  : " << mynewconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << mynewconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << mynewconfidence->GetExpectedCLb_b()  << endl;

// show canonical -2lnQ plots in a new canvas
// - The histogram of -2lnQ for background hypothesis (full)
// - The histogram of -2lnQ for signal and background hypothesis (dashed)
  TCanvas *c2 = new TCanvas("c2");
  myconfidence->Draw();
  
// clean up (except histograms and canvas)
  delete myconfidence;
//.........这里部分代码省略.........

void MakeADRecoParamEntry(AliRecoParam::EventSpecie_t defaultEventSpecie=AliRecoParam::kLowMult, const char *outputCDB = "local://$ALICE_ROOT/OCDB") {
//========================================================================
//
// Steering macro for AD reconstruction parameters
//
// Author: Michal Broz
//
//========================================================================

  const char* macroname = "MakeADRecoParam.C";

  // Activate CDB storage and load geometry from CDB
  AliCDBManager* cdb = AliCDBManager::Instance();
  cdb->SetDefaultStorage(outputCDB);
  cdb->SetRun(0);
  
  TObjArray *recoParamArray = new TObjArray();

  {
    AliADRecoParam * ADRecoParam = new AliADRecoParam;
    ADRecoParam->SetEventSpecie(AliRecoParam::kCosmic);
    ADRecoParam->SetStartClock(0);
    ADRecoParam->SetEndClock(20);
    ADRecoParam->SetNPreClocks(1);
    ADRecoParam->SetNPostClocks(10);
    ADRecoParam->SetTimeWindowBBALow(-2.5);
    ADRecoParam->SetTimeWindowBBAUp(2.5);
    ADRecoParam->SetTimeWindowBGALow(-4.0);
    ADRecoParam->SetTimeWindowBGAUp(4.0);
    ADRecoParam->SetTimeWindowBBCLow(-1.5);
    ADRecoParam->SetTimeWindowBBCUp(1.5);
    ADRecoParam->SetTimeWindowBGCLow(-2.0);
    ADRecoParam->SetTimeWindowBGCUp(2.0);
    ADRecoParam->SetAdcThresHold(5);
    ADRecoParam->SetMaxResid(1.5);
    ADRecoParam->SetResidRise(0.02);
    recoParamArray->AddLast(ADRecoParam);
  }
  {
    AliADRecoParam * ADRecoParam = new AliADRecoParam;
    ADRecoParam->SetStartClock(0);
    ADRecoParam->SetEndClock(20);
    ADRecoParam->SetNPreClocks(1);
    ADRecoParam->SetNPostClocks(10);
    ADRecoParam->SetTimeWindowBBALow(-2.5);
    ADRecoParam->SetTimeWindowBBAUp(2.5);
    ADRecoParam->SetTimeWindowBGALow(-4.0);
    ADRecoParam->SetTimeWindowBGAUp(4.0);
    ADRecoParam->SetTimeWindowBBCLow(-1.5);
    ADRecoParam->SetTimeWindowBBCUp(1.5);
    ADRecoParam->SetTimeWindowBGCLow(-2.0);
    ADRecoParam->SetTimeWindowBGCUp(2.0);
    ADRecoParam->SetAdcThresHold(5);
    ADRecoParam->SetMaxResid(1.5);
    ADRecoParam->SetResidRise(0.02);
    ADRecoParam->SetEventSpecie(AliRecoParam::kLowMult);
    recoParamArray->AddLast(ADRecoParam);
  }
  {
    AliADRecoParam * ADRecoParam = new AliADRecoParam;
    ADRecoParam->SetStartClock(9);
    ADRecoParam->SetEndClock(11);
    ADRecoParam->SetNPostClocks(6);
    ADRecoParam->SetTimeWindowBBALow(-2.5);
    ADRecoParam->SetTimeWindowBBAUp(2.5);
    ADRecoParam->SetTimeWindowBGALow(-4.0);
    ADRecoParam->SetTimeWindowBGAUp(4.0);
    ADRecoParam->SetTimeWindowBBCLow(-1.5);
    ADRecoParam->SetTimeWindowBBCUp(1.5);
    ADRecoParam->SetTimeWindowBGCLow(-2.0);
    ADRecoParam->SetTimeWindowBGCUp(2.0);
    ADRecoParam->SetAdcThresHold(5);
    ADRecoParam->SetMaxResid(1.5);
    ADRecoParam->SetResidRise(0.02);
    ADRecoParam->SetEventSpecie(AliRecoParam::kHighMult);
    recoParamArray->AddLast(ADRecoParam);
  }

  // Set the defaultEventSpecie
  Bool_t defaultIsSet = kFALSE;
  for(Int_t i =0; i < recoParamArray->GetEntriesFast(); i++) {
    AliDetectorRecoParam *param = (AliDetectorRecoParam *)recoParamArray->UncheckedAt(i);
    if (!param) continue;
    if (defaultEventSpecie & param->GetEventSpecie()) {
      param->SetAsDefault();
      defaultIsSet = kTRUE;
    }
  }

  if (!defaultIsSet) {
    Error(macroname,"The default reconstruction parameters are not set! Exiting...");
    return;
  }

  // save in CDB storage
  AliCDBMetaData *md= new AliCDBMetaData();
  md->SetResponsible("Michal Broz");
  md->SetComment("Reconstruction parameters for AD");
  md->SetAliRootVersion(gSystem->Getenv("ARVERSION"));
  md->SetBeamPeriod(0);
//.........这里部分代码省略.........

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

  AliCFTrackCutPid* cutPID = new AliCFTrackCutPid("cutPID","ESD_PID") ;
  int n_species = AliPID::kSPECIES ;
  Double_t* prior = new Double_t[n_species];
  
  prior[0] = 0.0244519 ;
  prior[1] = 0.0143988 ;
  prior[2] = 0.805747  ;
  prior[3] = 0.0928785 ;
  prior[4] = 0.0625243 ;
  
  cutPID->SetPriors(prior);
  cutPID->SetProbabilityCut(0.0);
  if (readTPCTracks) cutPID->SetDetectors("TPC");
  else               cutPID->SetDetectors("ALL");
  if (readAOD) cutPID->SetAODmode(kTRUE );
  else         cutPID->SetAODmode(kFALSE);
  switch(TMath::Abs(PDG)) {
  case 11   : cutPID->SetParticleType(AliPID::kElectron, kTRUE); break;
  case 13   : cutPID->SetParticleType(AliPID::kMuon    , kTRUE); break;
  case 211  : cutPID->SetParticleType(AliPID::kPion    , kTRUE); break;
  case 321  : cutPID->SetParticleType(AliPID::kKaon    , kTRUE); break;
  case 2212 : cutPID->SetParticleType(AliPID::kProton  , kTRUE); break;
  default   : printf("UNDEFINED PID\n"); break;
  }
  cutPID->SetQAOn(qaList);

  printf("CREATE EVENT LEVEL CUTS\n");
  TObjArray* evtList = new TObjArray(0) ;
//   evtList->AddLast(evtRecCuts);
  
  printf("CREATE MC KINE CUTS\n");
  TObjArray* mcList = new TObjArray(0) ;
  mcList->AddLast(mcKineCuts);
  mcList->AddLast(mcGenCuts);

  printf("CREATE ACCEPTANCE CUTS\n");
  TObjArray* accList = new TObjArray(0) ;
  accList->AddLast(mcAccCuts);

  printf("CREATE RECONSTRUCTION CUTS\n");
  TObjArray* recList = new TObjArray(0) ;
  recList->AddLast(recKineCuts);
  recList->AddLast(recQualityCuts);
  recList->AddLast(recIsPrimaryCuts);

  printf("CREATE PID CUTS\n");
  TObjArray* fPIDCutList = new TObjArray(0) ;
  fPIDCutList->AddLast(cutPID);

  //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
  printf("CREATE INTERFACE AND CUTS\n");
  AliCFManager* man = new AliCFManager() ;

  man->SetNStepEvent(1);
  man->SetEventCutsList(0,evtList);

  man->SetParticleContainer(container);
  man->SetParticleCutsList(0,mcList);
  man->SetParticleCutsList(1,accList);
  man->SetParticleCutsList(2,recList);
  man->SetParticleCutsList(3,fPIDCutList);


  //CREATE THE TASK
  printf("CREATE TASK\n");

void bfcread_tagsBranch( 
  const char *MainFile="/afs/rhic.bnl.gov/star/data/samples/gstar.tags.root",
  Int_t printEvent=1,
  const char *fname="qa_tags.out",
  Int_t fullPrint=0) 
{
  // start timer
  TStopwatch timer;
  timer.Start();

  cout << endl << endl;
  cout << " bfcread_tagsBranch.C: input file  = " << MainFile << endl;
  cout << " bfcread_tagsBranch.C: print event # " << printEvent << endl;
  cout << " bfcread_tagsBranch.C: output file = " << fname << endl;
  cout << " bfcread_tagsBranch.C: full printout = " << fullPrint << endl;
  cout << endl;

  ofstream fout(fname);

  fout << endl << endl;
  fout << " bfcread_tagsBranch.C: input file  = " << MainFile << endl;
  fout << " bfcread_tagsBranch.C: print evt#  = " << printEvent << endl;
  fout << " bfcread_tagsBranch.C: output file = " << fname << endl;
  fout << " bfcread_tagsBranch.C: full printout = " << fullPrint << endl;
  fout << endl;

  TFile *file = TFile::Open(MainFile);
  TTree *tree = (TTree*)file->Get("Tag");

  cout <<" read file: " << file->GetName() << endl << endl;

  Int_t nEntries = tree->GetEntries();
  cout << " Total # events  = " << nEntries << endl;

  TObjArray *leaves = tree->GetListOfLeaves();
  Int_t nLeaves = leaves->GetEntriesFast();

  cout << "  Total # leaves  = " << nLeaves << endl;

  TString *tName = new TString(" ");
  TNamed *tableName=0;
  TObjArray *tagTable = new TObjArray;
  Int_t tableCount = 0;
  Int_t *tableIndex = new Int_t[nLeaves];
  Int_t tagCount = 0;

  TBranch *branch=0;
  TLeaf *leaf=0;
  Int_t ndim =0;

//count number of tag tables encoded in the TTree branch names
  for (Int_t l=0;l<nLeaves;l++) {
    leaf = (TLeaf*)leaves->UncheckedAt(l);
    tagCount+=leaf->GetNdata();
    branch = leaf->GetBranch();
      cout << "leaf #  " << l << "  br name = " <<  
                                     branch->GetName() << endl;
    //new tag table name
    if ( strstr(branch->GetName(), tName->Data()) == 0 ) {
      tName = new TString(branch->GetName());
      tName->Resize(tName->Index("."));
      //the tableName is encoded in the branch Name before the "."
      tableName = new TNamed(tName->Data(),"Tag");
      tagTable->AddLast(tableName);
      tableCount++;
    }
    tableIndex[l]=tableCount-1;
  }

  cout << endl << "  Total num tables(branches),tags = " 
              << tableCount << "   " << tagCount << endl << endl;

  Int_t *countTagsTable = new Int_t[tableCount];
  Int_t *countLeavesTable = new Int_t[tableCount];
  Float_t *sumTagsLeaf = new Float_t[nLeaves];
  Int_t *countTagsLeaf = new Int_t[nLeaves];

// Now loop over leaves (to get values of tags)

   Int_t setBranch = -1;
   Int_t nowBranch = -1;

   for (Int_t l=0;l<nLeaves;l++) {

      leaf = (TLeaf*)leaves->UncheckedAt(l);
      branch = leaf->GetBranch();
      ndim = leaf->GetNdata();

      nowBranch = tableIndex[l];

      //cout << " nowbranch, setBranch = " << 
      //  nowBranch << ", "<< setBranch << endl;

      Float_t RtableIndex=tableIndex[l];

      if (nowBranch !=  setBranch){ 
          setBranch=nowBranch;
          cout << " QAInfo: branch ";
            cout.width(2);
            cout << tableIndex[l] << " = ";
//.........这里部分代码省略.........