PWG4/macros/AddTaskJets.C

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