[u/mrichter/AliRoot.git] / PWGJE / macros / AddTaskJets.C

AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask); // Common config
AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius = -1);

AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf,Float_t radius = -1,UInt_t filterMask = 0,Float_t ptTrackMin = 0.15,Int_t iBack = 0); // for the new AF
Int_t AddTaskJetsDelta(char *nonStdFile = "",UInt_t filterMask = 0,Bool_t kUseAODMC = kTRUE,UInt_t runFlag = 1|4|32|128|256);     
AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask = 0);

Float_t kPtTrackMin = 0.15;
Int_t kBackgroundMode = 0;


AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){
  // fills the standard "jets" branch in the AOD
  // need the ESDFilter to run before, to access the AODtracks
  // Tracks selected by the first Filter (1<<0)
  // needs to be adapted for different cuts
  
  // UA1 as standard chosen, since it is the most robust and simple JF
  // R = 0.4 suffficient to provide accurate jet axis for correlation studies
  // energy resolution suffers a little
  // Acceptance of jets not limited by the Jet Finder but should be done
  // by user to abs(eta) < 0.5 

  return AddTaskJets("AOD","UA1",0.4,filterMask);

}


Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){

  // Adds a whole set of jet finders  all to be written
  // to a delta AOD
  
  // this can in principle be done also with on the fly 
  // if we have an ESD input jet fidner task does automatically fetch the ouput aod

   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   if (!mgr) {
      ::Error("AddTaskJetsDelta", "No analysis manager to connect to.");
      return 0;
   }  
   
   // Check the analysis type using the event handlers connected to the analysis manager.
   //==============================================================================
   if (!mgr->GetInputEventHandler()) {
      ::Error("AddTaskJetsDelta", "This task requires an input event handler");
      return 0;
   }

  AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();
  if (!aodh) {
    ::Error("AddTaskJetsDelta", "This task needs an output event handler");
    return 0;
  }   

  TString outFile(nonStdFile);

  AliAnalysisTaskJets *jetana = 0;
  Int_t iCount = 0;

  // Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc
  const char *cJF[14]        = {"UA1","UA1","UA1", "CDF", "CDF","DA","DA","SISCONE","SISCONE","FASTJET","FASTJET","FASTKT","FASTKT","UA1LO"};
  const Float_t radius[14]   = {  0.4,  0.7,  1.0,  0.4,    0.7, 0.4, 0.7,      0.4,     0.7,      0.4,       0.7,     0.4,     0.7,    0.4};
  UInt_t  flag[14]           = {    6,    7,    7,    7,      7,   7,   7,        7,       7,        7,         7,       7,       7,      7};
  // flag[5] = 0; // set siscone to 0 for proof mode...
  // flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)
  // i.e. 7 all, 6 only MC2 and MC
  // this stay at three
  const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"};  

  
  for(int i = 0; i< 14;i++){
    if(!(runFlag&(1<<i)))continue;
    if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it
    for(int ib = 0;ib<4;ib++){      
      if(flag[i]&(1<<ib)){
	jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);
	if(jetana){
	  char *cRadius = "";
	  if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision
	  //	  jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets
	  if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());
	  iCount++;
	}
      }
    }
  }
    
  Printf("Added %d JetFinders",iCount);
}


AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf, Float_t radius,UInt_t filterMask,float ptTrackMin,int iBack)
{
  // Creates a jet finder task, configures it and adds it to the analysis manager.

  kPtTrackMin = ptTrackMin;
  kBackgroundMode = iBack;
  // Get the pointer to the existing analysis manager via the static access method.
  //==============================================================================
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  if (!mgr) {
    ::Error("AddTaskJets", "No analysis manager to connect to.");
    return NULL;
  }  
  
  // Check the analysis type using the event handlers connected to the analysis manager.
  //==============================================================================
  if (!mgr->GetInputEventHandler()) {
    ::Error("AddTaskJets", "This task requires an input event handler");
    return NULL;
   }
  
  AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();
  if (!aodh) {
    ::Error("AddTaskJets", "This task needs an output event handler");
    return NULL;
  }   
  

  // Create the task and configure it.
  //===========================================================================
  AliAnalysisTaskJets *jetana;
   AliJetReader *er = CreateJetReader(jr,filterMask);
   // Define jet header and jet finder
   AliJetFinder *jetFinder = CreateJetFinder(jf,radius);

   if (jetFinder){
       if (er) jetFinder->SetJetReader(er);
   }


   TString cAdd = "";
   cAdd += Form("%02d",(int)((radius+0.01)*10.));
   cAdd += Form("_B%d",(int)((kBackgroundMode)));
   cAdd += Form("_Filter%05d",filterMask);
   cAdd += Form("_Cut%05d",(int)((1000.*kPtTrackMin)));
   Printf("%s",cAdd.Data());
   jetana = new AliAnalysisTaskJets(Form("JetAnalysis%s_%s%s",jr,jf,cAdd.Data()));


   TString c_jr(jr);
   c_jr.ToLower();
   TString c_jf(jf);
   c_jf.ToLower();

   /*
   if(c_jf.CompareTo("ua1")==0 && TMath::Abs(radius-0.4) < 0.01 && c_jr.CompareTo("aod") == 0){
     // do nothing, this is the standard jet finder R = 0.4, UA1 on AOD
   }
   else{

     
 }
   */
   TString bName =  Form("jets%s_%s%s",jr,jf,cAdd.Data());
   jetana->SetNonStdBranch(bName.Data());
   Printf("Set jet branchname \"%s\"",bName.Data());


   AliAnalysisDataContainer *cout_jet = mgr->CreateContainer(
							     Form("jethist_%s_%s%s",
								  c_jr.Data(),
								  c_jf.Data(),
								  cAdd.Data()), 
							     TList::Class(),
							     AliAnalysisManager::kOutputContainer, 
							     Form("%s:PWGJE_jethist_%s_%s%s",AliAnalysisManager::GetCommonFileName(),
								  c_jr.Data(),
								  c_jf.Data(),
								  cAdd.Data()));

   // Connect jet finder to task.
   jetana->SetJetFinder(jetFinder);
   jetana->SetConfigFile("");
   jetana->SetDebugLevel(2);
   if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){
     //     jetana->SetFilterPt(20.);
   }


   mgr->AddTask(jetana);


   // Create ONLY the output containers for the data produced by the task.
   // Get and connect other common input/output containers via the manager as below
   //==============================================================================
   mgr->ConnectInput  (jetana, 0, mgr->GetCommonInputContainer());
// AOD output slot will be used in a different way in future
   mgr->ConnectOutput (jetana, 0, mgr->GetCommonOutputContainer());
   mgr->ConnectOutput (jetana, 1, cout_jet);
   
   return jetana;
}

AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius){
  AliJetFinder *jetFinder = 0;

  switch (jf) {
  case "CDF":
    AliCdfJetHeader *jh = new AliCdfJetHeader();
    jh->SetRadius(0.7);
    if(radius>0)jh->SetRadius(radius);    
    jh->SetAODwrite(kTRUE);
    jh->SetAODtracksWrite(kTRUE);
    //    jh->SetDebugCDF(1);
    jetFinder = new AliCdfJetFinder();
    if (jh) jetFinder->SetJetHeader(jh);
    break;

  case "DA":
    AliDAJetHeader *jh=new AliDAJetHeader();
    jh->SetComment("DA jet code with default parameters");
    jh->SelectJets(kTRUE);
//    jh->SetNeff(200);
//    jh->SetEtaEff(2.2);
    if(radius>0)jh->SetRadius(radius);
    jh->SetEtMin(5.);
    jetFinder = new AliDAJetFinder();
    if (jh) jetFinder->SetJetHeader(jh);
    break;

  case "FASTJET":
    // DEFAULT is ANTI KT
    AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
    jh->SetRparam(0.4); // setup parameters                                  
    Double_t rBkg = 0.6;
    if(radius>0)jh->SetRparam(radius);
    jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh
    /*
      $FASTJET/include/fastjet/JetDefinition.hh
      enum JetAlgorithm {kt_algorithm, cambridge_algorithm,
      antikt_algorithm, genkt_algorithm,
      ee_kt_algorithm, ee_genkt_algorithm, ...};

     */
    jetFinder = new AliFastJetFinder();
    jh->SetPtMin(1);
    jh->SetRparamBkg(rBkg);
    jh->SetGhostEtaMax(0.9+rBkg);
    if (jh) jetFinder->SetJetHeader(jh);
    break;

  case "FASTKT":
    AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
    jh->SetRparam(0.4); // setup parameters                                  
    Double_t rBkg = 0.6;
    if(radius>0)jh->SetRparam(radius);
    jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh
    jh->SetPtMin(1);
    jh->SetRparamBkg(rBkg);
    jh->SetGhostEtaMax(0.9+rBkg);
    jetFinder = new AliFastJetFinder();
    if (jh) jetFinder->SetJetHeader(jh);
    break;

  case "SISCONE":
    AliSISConeJetHeader * jh = new AliSISConeJetHeader();

    jh->SetJetEtaMax(1.5);
    jh->SetJetEtaMin(-1.5);

    //siscone parameters
    jh->SetConeRadius(0.4);                   // default cone radius
    if(radius>0)jh->SetConeRadius(radius);   // cone radius

    jh->SetOverlapThreshold(0.75);            // overlap parameter, between 0 and 1 excluded!! 0.75 value is advised
    jh->SetPtProtojetMin(0);                  // pt min of protojets
    jh->SetMinJetPt(5);                      // Ptmin of jets (GeV)

    //do you want to subtract BG (0 = no, 1 = yes)
    jh->SetBGMode(0); // if 1 set also the radius for the background determination..

    //for background
    jh->SetRapRange( -0.9, 0.9);              // rapidity range for subtracting background must be < ghostmaxrap-0.95*R
    jh->SetPhiRange(0 , 2*TMath::Pi());       // phi range for subtracting background
    
    //to determine jets area
    jh->SetBGAlgorithm(1);                    // algorithm for rho calculation : 0 = kT, 1 = Cambridge
    jh->SetGhostEtaMax(4);                    // eta max where ghosts are generated 
    jh->SetGhostArea(0.05);                   // area of a ghost 
    jh->SetMeanGhostKt(1e-100);               // average transverse momentum of the ghosts.
    jh->SetAreaTypeNumber(4);                 // from 1 to 5 : 1 = active_area, 2 = active_area_explicit_ghosts, 3 = one_ghost_passive_area, 4 = passive_area, 5 = voronoi_area
    jetFinder = new AliSISConeJetFinder();
    if (jh) jetFinder->SetJetHeader(jh);
    break;

  case "UA1":
    AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
    jh->SetComment("UA1 jet code with default parameters");
    jh->BackgMode(0);
    jh->SetRadius(0.4);
    if(radius>0)jh->SetRadius(radius);
    jh->SetEtSeed(4.);
    jh->SetNAcceptJets(6);
    jh->SetLegoNbinPhi(432);
    jh->SetLegoNbinEta(274);
    jh->SetLegoEtaMin(-2);
    jh->SetLegoEtaMax(+2);
    jh->SetMinJetEt(5.);
    jh->SetJetEtaMax(1.5);
    jh->SetJetEtaMin(-1.5);
    jh->BackgMode(kBackgroundMode);

    jetFinder = new AliUA1JetFinderV1();
    if (jh) jetFinder->SetJetHeader(jh);
    break;
  case "UA1LO":
    AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
    jh->SetComment("UA1 jet code with Lo Pt settings parameters");
    jh->BackgMode(0);
    jh->SetRadius(0.4);
    if(radius>0)jh->SetRadius(radius);
    jh->SetEtSeed(1.);
    jh->SetNAcceptJets(6);
    jh->SetLegoNbinPhi(432);
    jh->SetLegoNbinEta(274);
    jh->SetLegoEtaMin(-2);
    jh->SetLegoEtaMax(+2);
    jh->SetMinJetEt(1.);
    jh->SetJetEtaMax(1.5);
    jh->SetJetEtaMin(-1.5);
    jh->BackgMode(kBackgroundMode);

    jetFinder = new AliUA1JetFinderV1();
    if (jh) jetFinder->SetJetHeader(jh);
    break;


  case "UA1MC":
    AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
    jh->SetComment("UA1 jet code with default MC parameters");
    jh->BackgMode(0);
    jh->SetRadius(1.0);
    if(radius>0)jh->SetRadius(radius);
    jh->SetEtSeed(4.);
    jh->SetNAcceptJets(6);
    jh->SetLegoNbinPhi(432);
    jh->SetLegoNbinEta(274);
    jh->SetLegoEtaMin(-2);
    jh->SetLegoEtaMax(+2);
    jh->SetMinJetEt(5.);
    jh->SetJetEtaMax(1.5);
    jh->SetJetEtaMin(-1.5);
    jh->BackgMode(kBackgroundMode);
    jetFinder = new AliUA1JetFinderV1();
    if (jh) jetFinder->SetJetHeader(jh);
    break;
  default:
    Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");
    break;
  }

  return jetFinder;

}

AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask){
  AliJetReader *er = 0;
  Float_t ptCut  = kPtTrackMin ; // cut on track p_T
  switch (jr) {
  case "MC":
    AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
    jrh->SetComment("MC full Kinematics");
    jrh->SetFastSimTPC(kFALSE);
    jrh->SetFastSimEMCAL(kFALSE);
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9                                                                             
    // Define reader and set its header                                     
    er = new AliJetKineReader();
    er->SetReaderHeader(jrh);
    break;
  case "MC2":
    AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
    jrh->SetComment("MC full Kinematics spearate config charged only");
    jrh->SetFastSimTPC(kFALSE);
    jrh->SetFastSimEMCAL(kFALSE);
    jrh->SetChargedOnly(kTRUE);
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9                                                                             
    // Define reader and set its header                                     
    er = new AliJetKineReader();
    er->SetReaderHeader(jrh);
    break;
  case "ESD":
    AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();
    jrh->SetComment("Testing");
    jrh->SetFirstEvent(0.);
    jrh->SetLastEvent(1000);
    jrh->SetPtCut(ptCut);
    jrh->SetReadSignalOnly(kFALSE);
    // Define reader and set its header                                     
    er = new AliJetESDReader();
    er->SetReaderHeader(jrh);
    break;

  case "AOD":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD Reader");
    jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
    jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
    if(filterMask>0)jrh->SetTestFilterMask(filterMask); 
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODMC":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD MC Reader");
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
    jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODMC2":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD MC Reader");
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
    jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODMC2b":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD MC Reader");
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-0.9,0.9); // to take all MC particles default is 0.9
    jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODextra":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD Reader with extra branch");
    jrh->SetPtCut(ptCut); // set low p_T cut of to ptCut
    jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
    jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
    jrh->SetNonStdBranch("aodExtraTracks");
    if(filterMask>0)jrh->SetTestFilterMask(filterMask); 
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODextraonly":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD Reader with extra branch");
    jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
    jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
    jrh->SetReadBranches(AliJetAODReaderHeader::kReadNonStdBranch);
    jrh->SetNonStdBranch("aodExtraTracks");
    if(filterMask>0)jrh->SetTestFilterMask(filterMask); 
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODMCextra":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD MC Reader with extra branch");
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
    jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
    jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
    jrh->SetNonStdBranch("aodExtraMCparticles");
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  case "AODMC2extra":
    AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
    jrh->SetComment("AOD MC Reader with extra branch");
    jrh->SetPtCut(ptCut);
    jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
    jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
    jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
    jrh->SetNonStdBranch("aodExtraMCparticles");
    // Define reader and set its header
    er = new AliJetAODReader();
    er->SetReaderHeader(jrh);
    break;
  default:
    ::Error("AddTaskJets", "Wrong jet reader selected\n");
    return 0;
  }

  return er;

}
Commit	Line	Data
a1d44e31	1	AliJetReader CreateJetReader(Char_t jr,UInt_t filterMask); // Common config
	2	AliJetFinder CreateJetFinder(Char_t jf,Float_t radius = -1);
	3
	4	AliAnalysisTaskJets AddTaskJets(Char_t jr, Char_t *jf,Float_t radius = -1,UInt_t filterMask = 0,Float_t ptTrackMin = 0.15,Int_t iBack = 0); // for the new AF
	5	Int_t AddTaskJetsDelta(char *nonStdFile = "",UInt_t filterMask = 0,Bool_t kUseAODMC = kTRUE,UInt_t runFlag = 1\|4\|32\|128\|256);
	6	AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask = 0);
	7
	8	Float_t kPtTrackMin = 0.15;
	9	Int_t kBackgroundMode = 0;
	10
	11
	12	AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){
	13	// fills the standard "jets" branch in the AOD
	14	// need the ESDFilter to run before, to access the AODtracks
	15	// Tracks selected by the first Filter (1<<0)
	16	// needs to be adapted for different cuts
	17
	18	// UA1 as standard chosen, since it is the most robust and simple JF
	19	// R = 0.4 suffficient to provide accurate jet axis for correlation studies
	20	// energy resolution suffers a little
	21	// Acceptance of jets not limited by the Jet Finder but should be done
	22	// by user to abs(eta) < 0.5
	23
	24	return AddTaskJets("AOD","UA1",0.4,filterMask);
	25
	26	}
	27
	28
	29
	30	Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){
	31
	32	// Adds a whole set of jet finders all to be written
	33	// to a delta AOD
	34
	35	// this can in principle be done also with on the fly
	36	// if we have an ESD input jet fidner task does automatically fetch the ouput aod
	37
	38	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	39	if (!mgr) {
	40	::Error("AddTaskJetsDelta", "No analysis manager to connect to.");
	41	return 0;
	42	}
	43
	44	// Check the analysis type using the event handlers connected to the analysis manager.
	45	//==============================================================================
	46	if (!mgr->GetInputEventHandler()) {
	47	::Error("AddTaskJetsDelta", "This task requires an input event handler");
	48	return 0;
	49	}
	50
	51	AliAODHandler aodh = (AliAODHandler)mgr->GetOutputEventHandler();
	52	if (!aodh) {
	53	::Error("AddTaskJetsDelta", "This task needs an output event handler");
	54	return 0;
	55	}
	56
	57	TString outFile(nonStdFile);
	58
	59	AliAnalysisTaskJets *jetana = 0;
	60	Int_t iCount = 0;
	61
	62	// Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc
	63	const char *cJF[14] = {"UA1","UA1","UA1", "CDF", "CDF","DA","DA","SISCONE","SISCONE","FASTJET","FASTJET","FASTKT","FASTKT","UA1LO"};
	64	const Float_t radius[14] = { 0.4, 0.7, 1.0, 0.4, 0.7, 0.4, 0.7, 0.4, 0.7, 0.4, 0.7, 0.4, 0.7, 0.4};
65	UInt_t flag[14] = { 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7};
66	// flag[5] = 0; // set siscone to 0 for proof mode...
67	// flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)
68	// i.e. 7 all, 6 only MC2 and MC
69	// this stay at three
70	const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"};
71
72
73
74	for(int i = 0; i< 14;i++){
75	if(!(runFlag&(1<<i)))continue;
76	if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it
77	for(int ib = 0;ib<4;ib++){
78	if(flag[i]&(1<<ib)){
79	jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);
80	if(jetana){
81	char *cRadius = "";
82	if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision
83	// jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets
84	if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());
85	iCount++;
86	}
87	}
88	}
89	}
90
91	Printf("Added %d JetFinders",iCount);
92	}
93
94
95
96
97
98	AliAnalysisTaskJets AddTaskJets(Char_t jr, Char_t *jf, Float_t radius,UInt_t filterMask,float ptTrackMin,int iBack)
99	{
100	// Creates a jet finder task, configures it and adds it to the analysis manager.
101
102	kPtTrackMin = ptTrackMin;
103	kBackgroundMode = iBack;
104	// Get the pointer to the existing analysis manager via the static access method.
105	//==============================================================================
106	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
107	if (!mgr) {
108	::Error("AddTaskJets", "No analysis manager to connect to.");
109	return NULL;
110	}
111
112	// Check the analysis type using the event handlers connected to the analysis manager.
113	//==============================================================================
114	if (!mgr->GetInputEventHandler()) {
115	::Error("AddTaskJets", "This task requires an input event handler");
116	return NULL;
117	}
118
119	AliAODHandler aodh = (AliAODHandler)mgr->GetOutputEventHandler();
120	if (!aodh) {
121	::Error("AddTaskJets", "This task needs an output event handler");
122	return NULL;
123	}
124
125
126	// Create the task and configure it.
127	//===========================================================================
128	AliAnalysisTaskJets *jetana;
129	AliJetReader *er = CreateJetReader(jr,filterMask);
130	// Define jet header and jet finder
131	AliJetFinder *jetFinder = CreateJetFinder(jf,radius);
132
133	if (jetFinder){
134	if (er) jetFinder->SetJetReader(er);
135	}
136
137
138
139	TString cAdd = "";
140	cAdd += Form("%02d",(int)((radius+0.01)*10.));
141	cAdd += Form("_B%d",(int)((kBackgroundMode)));
142	cAdd += Form("_Filter%05d",filterMask);
143	cAdd += Form("_Cut%05d",(int)((1000.*kPtTrackMin)));
144	Printf("%s",cAdd.Data());
145	jetana = new AliAnalysisTaskJets(Form("JetAnalysis%s_%s%s",jr,jf,cAdd.Data()));
146
147
148	TString c_jr(jr);
149	c_jr.ToLower();
150	TString c_jf(jf);
151	c_jf.ToLower();
152
153	/*
154	if(c_jf.CompareTo("ua1")==0 && TMath::Abs(radius-0.4) < 0.01 && c_jr.CompareTo("aod") == 0){
155	// do nothing, this is the standard jet finder R = 0.4, UA1 on AOD
156	}
157	else{
158
159
160	}
161	*/
162	TString bName = Form("jets%s_%s%s",jr,jf,cAdd.Data());
163	jetana->SetNonStdBranch(bName.Data());
164	Printf("Set jet branchname \"%s\"",bName.Data());
165
166
167	AliAnalysisDataContainer *cout_jet = mgr->CreateContainer(
168	Form("jethist_%s_%s%s",
169	c_jr.Data(),
170	c_jf.Data(),
171	cAdd.Data()),
172	TList::Class(),
173	AliAnalysisManager::kOutputContainer,
5c4489e2	174	Form("%s:PWGJE_jethist_%s_%s%s",AliAnalysisManager::GetCommonFileName(),
a1d44e31	175	c_jr.Data(),
	176	c_jf.Data(),
	177	cAdd.Data()));
	178
	179	// Connect jet finder to task.
	180	jetana->SetJetFinder(jetFinder);
	181	jetana->SetConfigFile("");
	182	jetana->SetDebugLevel(2);
	183	if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){
	184	// jetana->SetFilterPt(20.);
	185	}
	186
	187
	188	mgr->AddTask(jetana);
	189
	190
	191	// Create ONLY the output containers for the data produced by the task.
	192	// Get and connect other common input/output containers via the manager as below
	193	//==============================================================================
	194	mgr->ConnectInput (jetana, 0, mgr->GetCommonInputContainer());
	195	// AOD output slot will be used in a different way in future
	196	mgr->ConnectOutput (jetana, 0, mgr->GetCommonOutputContainer());
	197	mgr->ConnectOutput (jetana, 1, cout_jet);
	198
	199	return jetana;
	200	}
	201
	202	AliJetFinder CreateJetFinder(Char_t jf,Float_t radius){
	203	AliJetFinder *jetFinder = 0;
	204
	205	switch (jf) {
	206	case "CDF":
	207	AliCdfJetHeader *jh = new AliCdfJetHeader();
	208	jh->SetRadius(0.7);
	209	if(radius>0)jh->SetRadius(radius);
	210	jh->SetAODwrite(kTRUE);
	211	jh->SetAODtracksWrite(kTRUE);
	212	// jh->SetDebugCDF(1);
	213	jetFinder = new AliCdfJetFinder();
	214	if (jh) jetFinder->SetJetHeader(jh);
	215	break;
	216
	217	case "DA":
	218	AliDAJetHeader *jh=new AliDAJetHeader();
	219	jh->SetComment("DA jet code with default parameters");
	220	jh->SelectJets(kTRUE);
	221	// jh->SetNeff(200);
	222	// jh->SetEtaEff(2.2);
	223	if(radius>0)jh->SetRadius(radius);
	224	jh->SetEtMin(5.);
	225	jetFinder = new AliDAJetFinder();
	226	if (jh) jetFinder->SetJetHeader(jh);
	227	break;
	228
	229	case "FASTJET":
	230	// DEFAULT is ANTI KT
	231	AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
	232	jh->SetRparam(0.4); // setup parameters
	233	Double_t rBkg = 0.6;
	234	if(radius>0)jh->SetRparam(radius);
	235	jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh
	236	/*
	237	$FASTJET/include/fastjet/JetDefinition.hh
	238	enum JetAlgorithm {kt_algorithm, cambridge_algorithm,
239	antikt_algorithm, genkt_algorithm,
240	ee_kt_algorithm, ee_genkt_algorithm, ...};
241
242	*/
243	jetFinder = new AliFastJetFinder();
244	jh->SetPtMin(1);
245	jh->SetRparamBkg(rBkg);
246	jh->SetGhostEtaMax(0.9+rBkg);
247	if (jh) jetFinder->SetJetHeader(jh);
248	break;
249
250	case "FASTKT":
251	AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
252	jh->SetRparam(0.4); // setup parameters
253	Double_t rBkg = 0.6;
254	if(radius>0)jh->SetRparam(radius);
255	jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh
256	jh->SetPtMin(1);
257	jh->SetRparamBkg(rBkg);
258	jh->SetGhostEtaMax(0.9+rBkg);
259	jetFinder = new AliFastJetFinder();
260	if (jh) jetFinder->SetJetHeader(jh);
261	break;
262
263	case "SISCONE":
264	AliSISConeJetHeader * jh = new AliSISConeJetHeader();
265
266	jh->SetJetEtaMax(1.5);
267	jh->SetJetEtaMin(-1.5);
268
269	//siscone parameters
270	jh->SetConeRadius(0.4); // default cone radius
271	if(radius>0)jh->SetConeRadius(radius); // cone radius
272
273	jh->SetOverlapThreshold(0.75); // overlap parameter, between 0 and 1 excluded!! 0.75 value is advised
274	jh->SetPtProtojetMin(0); // pt min of protojets
275	jh->SetMinJetPt(5); // Ptmin of jets (GeV)
276
277	//do you want to subtract BG (0 = no, 1 = yes)
278	jh->SetBGMode(0); // if 1 set also the radius for the background determination..
279
280	//for background
281	jh->SetRapRange( -0.9, 0.9); // rapidity range for subtracting background must be < ghostmaxrap-0.95*R
282	jh->SetPhiRange(0 , 2*TMath::Pi()); // phi range for subtracting background
283
284	//to determine jets area
285	jh->SetBGAlgorithm(1); // algorithm for rho calculation : 0 = kT, 1 = Cambridge
286	jh->SetGhostEtaMax(4); // eta max where ghosts are generated
287	jh->SetGhostArea(0.05); // area of a ghost
288	jh->SetMeanGhostKt(1e-100); // average transverse momentum of the ghosts.
289	jh->SetAreaTypeNumber(4); // from 1 to 5 : 1 = active_area, 2 = active_area_explicit_ghosts, 3 = one_ghost_passive_area, 4 = passive_area, 5 = voronoi_area
290	jetFinder = new AliSISConeJetFinder();
291	if (jh) jetFinder->SetJetHeader(jh);
292	break;
293
294	case "UA1":
295	AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
296	jh->SetComment("UA1 jet code with default parameters");
297	jh->BackgMode(0);
298	jh->SetRadius(0.4);
299	if(radius>0)jh->SetRadius(radius);
300	jh->SetEtSeed(4.);
301	jh->SetNAcceptJets(6);
302	jh->SetLegoNbinPhi(432);
303	jh->SetLegoNbinEta(274);
304	jh->SetLegoEtaMin(-2);
305	jh->SetLegoEtaMax(+2);
306	jh->SetMinJetEt(5.);
307	jh->SetJetEtaMax(1.5);
308	jh->SetJetEtaMin(-1.5);
309	jh->BackgMode(kBackgroundMode);
310
311	jetFinder = new AliUA1JetFinderV1();
312	if (jh) jetFinder->SetJetHeader(jh);
313	break;
314	case "UA1LO":
315	AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
316	jh->SetComment("UA1 jet code with Lo Pt settings parameters");
317	jh->BackgMode(0);
318	jh->SetRadius(0.4);
319	if(radius>0)jh->SetRadius(radius);
320	jh->SetEtSeed(1.);
321	jh->SetNAcceptJets(6);
322	jh->SetLegoNbinPhi(432);
323	jh->SetLegoNbinEta(274);
324	jh->SetLegoEtaMin(-2);
325	jh->SetLegoEtaMax(+2);
326	jh->SetMinJetEt(1.);
327	jh->SetJetEtaMax(1.5);
328	jh->SetJetEtaMin(-1.5);
329	jh->BackgMode(kBackgroundMode);
330
331	jetFinder = new AliUA1JetFinderV1();
332	if (jh) jetFinder->SetJetHeader(jh);
333	break;
334
335
336
337	case "UA1MC":
338	AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
339	jh->SetComment("UA1 jet code with default MC parameters");
340	jh->BackgMode(0);
341	jh->SetRadius(1.0);
342	if(radius>0)jh->SetRadius(radius);
343	jh->SetEtSeed(4.);
344	jh->SetNAcceptJets(6);
345	jh->SetLegoNbinPhi(432);
346	jh->SetLegoNbinEta(274);
347	jh->SetLegoEtaMin(-2);
348	jh->SetLegoEtaMax(+2);
349	jh->SetMinJetEt(5.);
350	jh->SetJetEtaMax(1.5);
351	jh->SetJetEtaMin(-1.5);
352	jh->BackgMode(kBackgroundMode);
353	jetFinder = new AliUA1JetFinderV1();
354	if (jh) jetFinder->SetJetHeader(jh);
355	break;
356	default:
357	Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");
358	break;
359	}
360
361	return jetFinder;
362
363	}
364
365	AliJetReader CreateJetReader(Char_t jr,UInt_t filterMask){
366	AliJetReader *er = 0;
367	Float_t ptCut = kPtTrackMin ; // cut on track p_T
368	switch (jr) {
369	case "MC":
370	AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
371	jrh->SetComment("MC full Kinematics");
372	jrh->SetFastSimTPC(kFALSE);
373	jrh->SetFastSimEMCAL(kFALSE);
374	jrh->SetPtCut(ptCut);
375	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9
376	// Define reader and set its header
377	er = new AliJetKineReader();
378	er->SetReaderHeader(jrh);
379	break;
380	case "MC2":
381	AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
382	jrh->SetComment("MC full Kinematics spearate config charged only");
383	jrh->SetFastSimTPC(kFALSE);
384	jrh->SetFastSimEMCAL(kFALSE);
385	jrh->SetChargedOnly(kTRUE);
386	jrh->SetPtCut(ptCut);
387	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9
388	// Define reader and set its header
389	er = new AliJetKineReader();
390	er->SetReaderHeader(jrh);
391	break;
392	case "ESD":
393	AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();
394	jrh->SetComment("Testing");
395	jrh->SetFirstEvent(0.);
396	jrh->SetLastEvent(1000);
397	jrh->SetPtCut(ptCut);
398	jrh->SetReadSignalOnly(kFALSE);
399	// Define reader and set its header
400	er = new AliJetESDReader();
401	er->SetReaderHeader(jrh);
402	break;
403
404	case "AOD":
405	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
406	jrh->SetComment("AOD Reader");
407	jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
408	jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
409	if(filterMask>0)jrh->SetTestFilterMask(filterMask);
410	// Define reader and set its header
411	er = new AliJetAODReader();
412	er->SetReaderHeader(jrh);
413	break;
414	case "AODMC":
415	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
416	jrh->SetComment("AOD MC Reader");
417	jrh->SetPtCut(ptCut);
418	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
419	jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
420	// Define reader and set its header
421	er = new AliJetAODReader();
422	er->SetReaderHeader(jrh);
423	break;
424	case "AODMC2":
425	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
426	jrh->SetComment("AOD MC Reader");
427	jrh->SetPtCut(ptCut);
428	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
429	jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
430	// Define reader and set its header
431	er = new AliJetAODReader();
432	er->SetReaderHeader(jrh);
433	break;
434	case "AODMC2b":
435	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
436	jrh->SetComment("AOD MC Reader");
437	jrh->SetPtCut(ptCut);
438	jrh->SetFiducialEta(-0.9,0.9); // to take all MC particles default is 0.9
439	jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
440	// Define reader and set its header
441	er = new AliJetAODReader();
442	er->SetReaderHeader(jrh);
443	break;
444	case "AODextra":
445	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
446	jrh->SetComment("AOD Reader with extra branch");
447	jrh->SetPtCut(ptCut); // set low p_T cut of to ptCut
448	jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
449	jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch \| AliJetAODReaderHeader::kReadNonStdBranch);
450	jrh->SetNonStdBranch("aodExtraTracks");
451	if(filterMask>0)jrh->SetTestFilterMask(filterMask);
452	// Define reader and set its header
453	er = new AliJetAODReader();
454	er->SetReaderHeader(jrh);
455	break;
456	case "AODextraonly":
457	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
458	jrh->SetComment("AOD Reader with extra branch");
459	jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
460	jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
461	jrh->SetReadBranches(AliJetAODReaderHeader::kReadNonStdBranch);
462	jrh->SetNonStdBranch("aodExtraTracks");
463	if(filterMask>0)jrh->SetTestFilterMask(filterMask);
464	// Define reader and set its header
465	er = new AliJetAODReader();
466	er->SetReaderHeader(jrh);
467	break;
468	case "AODMCextra":
469	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
470	jrh->SetComment("AOD MC Reader with extra branch");
471	jrh->SetPtCut(ptCut);
472	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
473	jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
474	jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch \| AliJetAODReaderHeader::kReadNonStdBranch);
475	jrh->SetNonStdBranch("aodExtraMCparticles");
476	// Define reader and set its header
477	er = new AliJetAODReader();
478	er->SetReaderHeader(jrh);
479	break;
480	case "AODMC2extra":
481	AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
482	jrh->SetComment("AOD MC Reader with extra branch");
483	jrh->SetPtCut(ptCut);
484	jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
485	jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
486	jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch \| AliJetAODReaderHeader::kReadNonStdBranch);
487	jrh->SetNonStdBranch("aodExtraMCparticles");
488	// Define reader and set its header
489	er = new AliJetAODReader();
490	er->SetReaderHeader(jrh);
491	break;
492	default:
493	::Error("AddTaskJets", "Wrong jet reader selected\n");
494	return 0;
495	}
496
497	return er;
498
499	}