PWG4/macros/AddTaskJets.C

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