PWGJE/macros/AddTaskJetsFinder.C

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