PWGGA/EMCALTasks/macros/AddTaskESDFilterEMCALEventSelect.C

   1 // $Id$
   2
   3 Bool_t AddTrackCutsLHC10h(AliAnalysisTaskESDfilter* esdFilter);
   4 Bool_t AddTrackCutsLHC11h(AliAnalysisTaskESDfilter* esdFilter);
   5 Bool_t enableTPCOnlyAODTracksLocalFlag=kFALSE;
   6
   7 AliAnalysisTaskESDfilter *AddTaskESDFilterEMCALEventSelect(Float_t energyCut = 10,     // EMCAL
   8                                                            Int_t   ncellsCut = 2,      // EMCAL
   9                                                            Int_t   runNumber = 170000, // EMCAL
  10                                                            Bool_t  useKineFilter=kTRUE,
  11                                                            Int_t   tofTimeZeroType=AliESDpid::kTOF_T0,
  12                                                            Bool_t  enableTPCOnlyAODTracks=kFALSE,
  13                                                            Bool_t  disableCascades=kFALSE,
  14                                                            Bool_t  disableKinks=kFALSE,
  15                                                            Int_t   runFlag = 1100)
  16 {
  17   // Creates a filter task and adds it to the analysis manager.
  18   // Get the pointer to the existing analysis manager via the static access method.
  19   //==============================================================================
  20   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  21   if (!mgr) {
  22     ::Error("AddTaskESDFilter", "No analysis manager to connect to.");
  23     return NULL;
  24   }
  25
  26   // This task requires an ESD input handler and an AOD output handler.
  27   // Check this using the analysis manager.
  28   //===============================================================================
  29   TString type = mgr->GetInputEventHandler()->GetDataType();
  30   if (!type.Contains("ESD")) {
  31     ::Error("AddTaskESDFilter", "ESD filtering task needs the manager to have an ESD input handler.");
  32     return NULL;
  33   }
  34   // Check if AOD output handler exist.
  35   AliAODHandler *aod_h = (AliAODHandler*)mgr->GetOutputEventHandler();
  36   if (!aod_h) {
  37     ::Error("AddTaskESDFilter", "ESD filtering task needs the manager to have an AOD output handler.");
  38     return NULL;
  39   }
  40   // Check if MC handler is connected in case kine filter requested
  41   AliMCEventHandler *mcH = (AliMCEventHandler*)mgr->GetMCtruthEventHandler();
  42   if (!mcH && useKineFilter) {
  43     ::Error("AddTaskESDFilter", "No MC handler connected while kine filtering requested");
  44     return NULL;
  45   }
  46
  47   // Create the task, add it to the manager and configure it.
  48   //===========================================================================
  49   AliAnalysisTaskESDfilterEMCALEventSelect *esdfilter = new AliAnalysisTaskESDfilterEMCALEventSelect("ESD Filter : EMCAL event select");
  50
  51   esdfilter->DisableZDC();
  52   esdfilter->DisablePmdClusters();
  53
  54   // EMCAL settings
  55   esdfilter->SetEnergyCut(energyCut);
  56   esdfilter->SetNcellsCut(ncellsCut);
  57
  58   AliEMCALRecoUtils * reco = esdfilter->GetRecoUtils();
  59   reco->SwitchOnRejectExoticCluster();
  60
  61   // Pass the bad channels, need to access run number
  62   TString fileName="$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root";
  63   AliOADBContainer *contBC=new AliOADBContainer("");
  64   contBC->InitFromFile((char*)fileName.Data(),"AliEMCALBadChannels");
  65   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runNumber);
  66   if(arrayBC){
  67     TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject("pass1");
  68     if(arrayBCpass){
  69
  70       reco->SwitchOnBadChannelsRemoval();
  71       printf("*** EMCAL RecoUtils : REMOVE bad cells \n");
  72
  73       for (Int_t i=0; i<10; ++i) {
  74         TH2I *hbm = reco->GetEMCALChannelStatusMap(i);
  75         if (hbm)
  76           delete hbm;
  77         hbm=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
  78
  79         if (!hbm) {
  80           AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
  81           continue;
  82         }
  83
  84         hbm->SetDirectory(0);
  85         reco->SetEMCALChannelStatusMap(i,hbm);
  86       }
  87     } else printf("AliEMCALRecoUtils ---Do NOT remove bad channels 1\n");
  88   }  else  printf("AliEMCALRecoUtils ---Do NOT remove bad channels 2\n");
  89
  90   // From here keep sync with $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
  91
  92   esdfilter->SetTimeZeroType(tofTimeZeroType);
  93   if  (disableCascades) esdfilter->DisableCascades();
  94   if  (disableKinks)    esdfilter->DisableKinks();
  95
  96   mgr->AddTask(esdfilter);
  97
  98   // Filtering of MC particles (decays conversions etc)
  99   // this task has to go AFTER all other filter tasks
 100   // since it fills the AODMC array with all
 101   // selected MC Particles, only this way we have the
 102   // AODMCparticle information available for following tasks
 103   AliAnalysisTaskMCParticleFilter *kinefilter = 0;
 104   if (useKineFilter) {
 105     kinefilter = new AliAnalysisTaskMCParticleFilter("Particle Kine Filter");
 106     mgr->AddTask(kinefilter);
 107   }
 108
 109   enableTPCOnlyAODTracksLocalFlag = enableTPCOnlyAODTracks;
 110   if((runFlag/100)==10){
 111     AddTrackCutsLHC10h(esdfilter);
 112   }
 113   else {
 114     // default 11h
 115     AddTrackCutsLHC11h(esdfilter);
 116   }
 117
 118   // Filter with cuts on V0s
 119   AliESDv0Cuts*   esdV0Cuts = new AliESDv0Cuts("Standard V0 Cuts pp", "ESD V0 Cuts");
 120   esdV0Cuts->SetMinRadius(0.2);
 121   esdV0Cuts->SetMaxRadius(200);
 122   esdV0Cuts->SetMinDcaPosToVertex(0.05);
 123   esdV0Cuts->SetMinDcaNegToVertex(0.05);
 124   esdV0Cuts->SetMaxDcaV0Daughters(1.5);
 125   esdV0Cuts->SetMinCosinePointingAngle(0.99);
 126   AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter");
 127   v0Filter->AddCuts(esdV0Cuts);
 128
 129   esdfilter->SetV0Filter(v0Filter);
 130
 131   // Create ONLY the output containers for the data produced by the task.
 132   // Get and connect other common input/output containers via the manager as below
 133   //==============================================================================
 134   mgr->ConnectInput  (esdfilter,  0, mgr->GetCommonInputContainer());
 135   mgr->ConnectOutput (esdfilter,  0, mgr->GetCommonOutputContainer());
 136   if (useKineFilter) {
 137     mgr->ConnectInput  (kinefilter,  0, mgr->GetCommonInputContainer());
 138     mgr->ConnectOutput (kinefilter,  0, mgr->GetCommonOutputContainer());
 139     AliAnalysisDataContainer *coutputEx = mgr->CreateContainer("cFilterList", TList::Class(),
 140                                                                AliAnalysisManager::kOutputContainer,"pyxsec_hists.root");
 141     mgr->ConnectOutput (kinefilter,  1,coutputEx);
 142   }
 143   return esdfilter;
 144 }
 145
 146 Bool_t AddTrackCutsLHC10h(AliAnalysisTaskESDfilter* esdfilter)
 147 {
 148   Printf("%s%d: Creating Track Cuts for LH10h",(char*)__FILE__,__LINE__);
 149
 150   // Cuts on primary tracks
 151   AliESDtrackCuts* esdTrackCutsL = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 152
 153   // ITS stand-alone tracks
 154   AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("ITS stand-alone Track Cuts", "ESD Track Cuts");
 155   esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE);
 156
 157   // Pixel OR necessary for the electrons
 158   AliESDtrackCuts *itsStrong = new AliESDtrackCuts("ITSorSPD", "pixel requirement for ITS");
 159   itsStrong->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 160
 161
 162   // PID for the electrons
 163   AliESDpidCuts *electronID = new AliESDpidCuts("Electrons", "Electron PID cuts");
 164   electronID->SetTPCnSigmaCut(AliPID::kElectron, 3.);
 165
 166   // tighter cuts on primary particles for high pT tracks
 167   // take the standard cuts, which include already
 168   // ITSrefit and use only primaries...
 169
 170   // ITS cuts for new jet analysis
 171   //  gROOT->LoadMacro("$ALICE_ROOT/PWGJE/macros/CreateTrackCutsPWGJE.C");
 172   //  AliESDtrackCuts* esdTrackCutsHG0 = CreateTrackCutsPWGJE(10001006);
 173
 174   AliESDtrackCuts *jetCuts1006 = new AliESDtrackCuts("AliESDtrackCuts");
 175
 176   TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
 177   jetCuts1006->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
 178   jetCuts1006->SetMinNClustersTPC(70);
 179   jetCuts1006->SetMaxChi2PerClusterTPC(4);
 180   jetCuts1006->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
 181   jetCuts1006->SetAcceptKinkDaughters(kFALSE);
 182   jetCuts1006->SetRequireTPCRefit(kTRUE);
 183   jetCuts1006->SetMaxFractionSharedTPCClusters(0.4);
 184   // ITS
 185   jetCuts1006->SetRequireITSRefit(kTRUE);
 186   //accept secondaries
 187   jetCuts1006->SetMaxDCAToVertexXY(2.4);
 188   jetCuts1006->SetMaxDCAToVertexZ(3.2);
 189   jetCuts1006->SetDCAToVertex2D(kTRUE);
 190   //reject fakes
 191   jetCuts1006->SetMaxChi2PerClusterITS(36);
 192   jetCuts1006->SetMaxChi2TPCConstrainedGlobal(36);
 193
 194   jetCuts1006->SetRequireSigmaToVertex(kFALSE);
 195
 196   jetCuts1006->SetEtaRange(-0.9,0.9);
 197   jetCuts1006->SetPtRange(0.15, 1E+15.);
 198
 199   AliESDtrackCuts* esdTrackCutsHG0 = jetCuts1006->Clone("JetCuts10001006");
 200   esdTrackCutsHG0->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 201
 202   // throw out tracks with too low number of clusters in
 203   // the first pass (be consistent with TPC only tracks)
 204   // N.B. the number off crossed rows still acts on the tracks after
 205   // all iterations if we require tpc standalone, number of clusters
 206   // and chi2 TPC cuts act on track after the first iteration
 207   //   esdTrackCutsH0->SetRequireTPCStandAlone(kTRUE);
 208   //   esdTrackCutsH0->SetMinNClustersTPC(80); // <--- first pass
 209
 210   // the complement to the one with SPD requirement
 211   //  AliESDtrackCuts* esdTrackCutsHG1 = CreateTrackCutsPWGJE(10011006);
 212   AliESDtrackCuts* esdTrackCutsHG1 = jetCuts1006->Clone("JetCuts10011006");
 213   esdTrackCutsHG1->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
 214
 215   // the tracks that must not be taken pass this cut and
 216   // non HGC1 and HG
 217   //  AliESDtrackCuts* esdTrackCutsHG2 = CreateTrackCutsPWGJE(10021006);
 218   AliESDtrackCuts* esdTrackCutsHG2 = jetCuts1006->Clone("JetCuts10021006");
 219   esdTrackCutsHG2->SetMaxChi2PerClusterITS(1E10);
 220
 221   // standard cuts also used in R_AA analysis
 222   //   "Global track RAA analysis QM2011 + Chi2ITS<36";
 223   //  AliESDtrackCuts* esdTrackCutsH2 = CreateTrackCutsPWGJE(1000);
 224   AliESDtrackCuts* esdTrackCutsH2 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
 225   esdTrackCutsH2->SetMinNCrossedRowsTPC(120);
 226   esdTrackCutsH2->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 227   esdTrackCutsH2->SetMaxChi2PerClusterITS(36);
 228   esdTrackCutsH2->SetMaxFractionSharedTPCClusters(0.4);
 229   esdTrackCutsH2->SetMaxChi2TPCConstrainedGlobal(36);
 230
 231   esdTrackCutsH2->SetEtaRange(-0.9,0.9);
 232   esdTrackCutsH2->SetPtRange(0.15, 1e10);
 233
 234   //  AliESDtrackCuts* esdTrackCutsGCOnly = CreateTrackCutsPWGJE(10041006);
 235   AliESDtrackCuts* esdTrackCutsGCOnly = jetCuts1006->Clone("JetCuts10041006");
 236   esdTrackCutsGCOnly->SetRequireITSRefit(kFALSE);
 237
 238   // TPC only tracks
 239   AliESDtrackCuts* esdTrackCutsTPCCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 240   esdTrackCutsTPCCOnly->SetMinNClustersTPC(70);
 241
 242   // Compose the filter
 243   AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");
 244   // 1, 1<<0
 245   trackFilter->AddCuts(esdTrackCutsL);
 246   // 2 1<<1
 247   trackFilter->AddCuts(esdTrackCutsITSsa);
 248   // 4 1<<2
 249   trackFilter->AddCuts(itsStrong);
 250   itsStrong->SetFilterMask(1);        // AND with Standard track cuts
 251   // 8 1<<3
 252   trackFilter->AddCuts(electronID);
 253   electronID->SetFilterMask(4);       // AND with Pixel Cuts
 254   // 16 1<<4
 255   trackFilter->AddCuts(esdTrackCutsHG0);
 256   // 32 1<<5
 257   trackFilter->AddCuts(esdTrackCutsHG1);
 258   // 64 1<<6
 259   trackFilter->AddCuts(esdTrackCutsHG2);
 260   // 128 1<<7
 261   trackFilter->AddCuts(esdTrackCutsTPCCOnly); // add QM TPC only track cuts
 262   if(enableTPCOnlyAODTracksLocalFlag)esdfilter->SetTPCOnlyFilterMask(128);
 263   // 256 1<<8
 264   trackFilter->AddCuts(esdTrackCutsGCOnly);
 265   // 512 1<<9
 266   AliESDtrackCuts* esdTrackCutsHG1_tmp = new AliESDtrackCuts(*esdTrackCutsHG1); // avoid double delete
 267   trackFilter->AddCuts(esdTrackCutsHG1_tmp); // add once more for tpc only tracks
 268   // 1024 1<<10
 269   trackFilter->AddCuts(esdTrackCutsH2); // add r_aa cuts
 270
 271
 272
 273   esdfilter->SetGlobalConstrainedFilterMask(1<<8|1<<9); // these tracks are written out as global constrained tracks
 274   esdfilter->SetHybridFilterMaskGlobalConstrainedGlobal((1<<4)); // these normal global tracks will be marked as hybrid
 275   esdfilter->SetWriteHybridGlobalConstrainedOnly(kTRUE); // write only the complement
 276   //     esdfilter->SetTPCConstrainedFilterMask(1<<11); // these tracks are written out as tpc constrained tracks
 277
 278   esdfilter->SetTrackFilter(trackFilter);
 279   return kTRUE;
 280 }
 281
 282 Bool_t AddTrackCutsLHC11h(AliAnalysisTaskESDfilter* esdfilter)
 283 {
 284   Printf("%s%d: Creating Track Cuts LHC11h",(char*)__FILE__,__LINE__);
 285
 286   // Cuts on primary tracks
 287   AliESDtrackCuts* esdTrackCutsL = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 288
 289   // ITS stand-alone tracks
 290   AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("ITS stand-alone Track Cuts", "ESD Track Cuts");
 291   esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE);
 292
 293   // Pixel OR necessary for the electrons
 294   AliESDtrackCuts *itsStrong = new AliESDtrackCuts("ITSorSPD", "pixel requirement for ITS");
 295   itsStrong->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 296
 297
 298   // PID for the electrons
 299   AliESDpidCuts *electronID = new AliESDpidCuts("Electrons", "Electron PID cuts");
 300   electronID->SetTPCnSigmaCut(AliPID::kElectron, 3.);
 301
 302   // standard cuts with very loose DCA
 303   AliESDtrackCuts* esdTrackCutsH = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 304   esdTrackCutsH->SetMaxDCAToVertexXY(2.4);
 305   esdTrackCutsH->SetMaxDCAToVertexZ(3.2);
 306   esdTrackCutsH->SetDCAToVertex2D(kTRUE);
 307
 308   // standard cuts with tight DCA cut
 309   AliESDtrackCuts* esdTrackCutsH2 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 310
 311   // standard cuts with tight DCA but with requiring the first SDD cluster instead of an SPD cluster
 312   // tracks selected by this cut are exclusive to those selected by the previous cut
 313   AliESDtrackCuts* esdTrackCutsH3 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 314   esdTrackCutsH3->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
 315   esdTrackCutsH3->SetClusterRequirementITS(AliESDtrackCuts::kSDD, AliESDtrackCuts::kFirst);
 316
 317   // TPC only tracks: Optionally enable the writing of TPConly information
 318   // constrained to SPD vertex in the filter below
 319   AliESDtrackCuts* esdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 320   // The following line is needed for 2010 PbPb reprocessing and pp, but not for 2011 PbPb
 321   //esdTrackCutsTPCOnly->SetMinNClustersTPC(70);
 322
 323   // Extra cuts for hybrids
 324   // first the global tracks we want to take
 325   AliESDtrackCuts* esdTrackCutsHTG = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 326   esdTrackCutsHTG->SetName("Global Hybrid tracks, loose DCA");
 327   esdTrackCutsHTG->SetMaxDCAToVertexXY(2.4);
 328   esdTrackCutsHTG->SetMaxDCAToVertexZ(3.2);
 329   esdTrackCutsHTG->SetDCAToVertex2D(kTRUE);
 330   esdTrackCutsHTG->SetMaxChi2TPCConstrainedGlobal(36);
 331
 332   // Than the complementary tracks which will be stored as global
 333   // constraint, complement is done in the ESDFilter task
 334   AliESDtrackCuts* esdTrackCutsHTGC = new AliESDtrackCuts(*esdTrackCutsHTG);
 335   esdTrackCutsHTGC->SetName("Global Constraint Hybrid tracks, loose DCA no it requirement");
 336   esdTrackCutsHTGC->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
 337   esdTrackCutsHTGC->SetRequireITSRefit(kFALSE);
 338
 339   // standard cuts with tight DCA cut, using cluster cut instead of crossed rows (a la 2010 default)
 340   AliESDtrackCuts* esdTrackCutsH2Cluster = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE, 0);
 341
 342   // Compose the filter
 343   AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");
 344   // 1, 1<<0
 345   trackFilter->AddCuts(esdTrackCutsL);
 346   // 2, 1<<1
 347   trackFilter->AddCuts(esdTrackCutsITSsa);
 348   // 4, 1<<2
 349   trackFilter->AddCuts(itsStrong);
 350   itsStrong->SetFilterMask(1);        // AND with Standard track cuts
 351   // 8, 1<<3
 352   trackFilter->AddCuts(electronID);
 353   electronID->SetFilterMask(4);       // AND with Pixel Cuts
 354   // 16, 1<<4
 355   trackFilter->AddCuts(esdTrackCutsH);
 356   // 32, 1<<5
 357   trackFilter->AddCuts(esdTrackCutsH2);
 358   // 64, 1<<6
 359   trackFilter->AddCuts(esdTrackCutsH3);
 360   // 128 , 1 << 7
 361   trackFilter->AddCuts(esdTrackCutsTPCOnly);
 362   if(enableTPCOnlyAODTracksLocalFlag)esdfilter->SetTPCOnlyFilterMask(128);
 363   // 256, 1 << 8 Global Hybrids
 364   trackFilter->AddCuts(esdTrackCutsHTG);
 365   esdfilter->SetHybridFilterMaskGlobalConstrainedGlobal((1<<8)); // these normal global tracks will be marked as hybrid
 366   // 512, 1<< 9 GlobalConstraint Hybrids
 367   trackFilter->AddCuts(esdTrackCutsHTGC);
 368   esdfilter->SetGlobalConstrainedFilterMask(1<<9); // these tracks are written out as global constrained tracks
 369   esdfilter->SetWriteHybridGlobalConstrainedOnly(kTRUE); // write only the complement
 370   // 1024, 1<< 10
 371   trackFilter->AddCuts(esdTrackCutsH2Cluster);
 372   esdfilter->SetTrackFilter(trackFilter);
 373
 374   return kTRUE;
 375 }