PWG/EMCAL/AliEsdSkimTask.cxx

   1 // $Id$
   2 //
   3 // Task to skim ESD files.
   4 //
   5 // Author: C.Loizides
   6
   7 #include "AliEsdSkimTask.h"
   8 #include <TClonesArray.h>
   9 #include <TFile.h>
  10 #include <TTree.h>
  11 #include "AliAnalysisManager.h"
  12 #include "AliCentrality.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDFMD.h"
  15 #include "AliESDMuonTrack.h"
  16 #include "AliESDtrackCuts.h"
  17 #include "AliESDtrackCuts.h"
  18 #include "AliEsdTrackExt.h"
  19 #include "AliEventplane.h"
  20 #include "AliInputEventHandler.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliPicoTrack.h"
  23 #include "AliVEventHandler.h"
  24
  25 //_________________________________________________________________________________________________
  26 AliEsdSkimTask::AliEsdSkimTask(const char *opt) :
  27   AliAnalysisTaskSE(opt), fEvent(0), fTree(0), fCuts(0),
  28   fDoZDC(1), fDoV0(1), fDoT0(1), fDoTPCv(1), fDoSPDv(1), fDoPriv(1),
  29   fDoEmCs(1), fDoPCs(0), fDoEmT(1), fDoPT(0), fDoTracks(0), fDoFmd(1),
  30   fDoMult(1), fDoTof(0), fDoPileup(1), fDoClus(0), fDoMuonTracks(0), fEmcNames(""),
  31   fDoMiniTracks(0), fTracks("Tracks"), fPhosClusOnly(0), fEmcalClusOnly(0),
  32   fDoSaveBytes(0), fDoCent(1), fDoRP(1), fRemoveCP(0), fResetCov(1),
  33   fDoPicoTracks(0), fCheckCond(0)
  34 {
  35   // Constructor.
  36
  37   if (!opt)
  38     return;
  39
  40   fBranchNames = "ESD:AliESDHeader.,AliESDRun.";
  41
  42   DefineOutput(1, TTree::Class());
  43 }
  44
  45 //_________________________________________________________________________________________________
  46 void AliEsdSkimTask::UserExec(Option_t */*opt*/)
  47 {
  48   // Process event.
  49
  50   AliESDEvent *esdin = dynamic_cast<AliESDEvent*>(InputEvent());
  51   if (!esdin)
  52     return;
  53
  54   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  55
  56   fEvent->ResetStdContent();
  57
  58   TList* objsin  = esdin->GetList();
  59   TList* objsout = fEvent->GetList();
  60
  61   if (fCheckCond!=0) {
  62     Int_t mult = AliESDtrackCuts::GetReferenceMultiplicity(esdin, AliESDtrackCuts::kTrackletsITSTPC, 0.8);
  63     if (mult<10)
  64       return;
  65   }
  66
  67   AliESDHeader *header = dynamic_cast<AliESDHeader*>(objsout->FindObject("AliESDHeader"));
  68   if (header) {
  69     am->LoadBranch("AliESDHeader.");
  70     if (!fDoSaveBytes) {
  71       *header = *esdin->GetHeader();
  72     } else {
  73       AliESDHeader *hin = esdin->GetHeader();
  74       header->SetTriggerMask(hin->GetTriggerMask());
  75       header->SetOrbitNumber(hin->GetOrbitNumber());
  76       header->SetTimeStamp(hin->GetTimeStamp());
  77       header->SetEventType(hin->GetEventType());
  78       header->SetEventSpecie(hin->GetEventSpecie());
  79       header->SetEventNumberInFile(hin->GetEventNumberInFile());
  80       header->SetBunchCrossNumber(hin->GetBunchCrossNumber());
  81       header->SetPeriodNumber(hin->GetPeriodNumber());
  82       header->SetTriggerCluster(hin->GetTriggerCluster());
  83       header->SetL0TriggerInputs(hin->GetL0TriggerInputs());
  84       header->SetL1TriggerInputs(hin->GetL1TriggerInputs());
  85       header->SetL2TriggerInputs(hin->GetL2TriggerInputs());
  86       for (Int_t i=0;i<24;++i) {
  87         const char *name = hin->GetTriggerInputName(i,0);
  88         if (name)
  89           header->SetActiveTriggerInputs(name,i);
  90       }
  91       for (Int_t i=0;i<24;++i) {
  92         const char *name = hin->GetTriggerInputName(i,1);
  93         if (name)
  94           header->SetActiveTriggerInputs(name,i+24);
  95       }
  96       for (Int_t i=0;i<12;++i) {
  97         const char *name = hin->GetTriggerInputName(i,2);
  98         if (name)
  99           header->SetActiveTriggerInputs(name,i+48);
 100       }
 101     }
 102     header->SetUniqueID(((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected());
 103     if (am->GetInputEventHandler()) {
 104       TTree *tree = am->GetInputEventHandler()->GetTree();
 105       if (tree) {
 106         TFile *cfile = tree->GetCurrentFile();
 107         if (cfile) {
 108           TString cname(cfile->GetName());
 109           header->SetTitle(cname);
 110         }
 111       }
 112     }
 113   }
 114   AliESDRun *run = dynamic_cast<AliESDRun*>(objsout->FindObject("AliESDRun"));
 115   if (run) {
 116     am->LoadBranch("AliESDRun.");
 117     *run = *esdin->GetESDRun();
 118   }
 119   AliCentrality *cent = dynamic_cast<AliCentrality*>(objsout->FindObject("Centrality"));
 120   if (cent) {
 121     cent->Reset();
 122     AliCentrality *centin = esdin->GetCentrality();
 123     cent->SetQuality(centin->GetQuality());
 124     cent->SetCentralityV0M(centin->GetCentralityPercentileUnchecked("V0M"));
 125     cent->SetCentralityFMD(centin->GetCentralityPercentileUnchecked("FMD"));
 126     cent->SetCentralityTRK(centin->GetCentralityPercentileUnchecked("TRK"));
 127     cent->SetCentralityTKL(centin->GetCentralityPercentileUnchecked("TKL"));
 128     cent->SetCentralityCL0(centin->GetCentralityPercentileUnchecked("CL0"));
 129     cent->SetCentralityCL1(centin->GetCentralityPercentileUnchecked("CL1"));
 130     cent->SetCentralityV0MvsFMD(centin->GetCentralityPercentileUnchecked("V0MvsFMD"));
 131     cent->SetCentralityTKLvsV0M(centin->GetCentralityPercentileUnchecked("TKLvsV0M"));
 132     cent->SetCentralityZEMvsZDC(centin->GetCentralityPercentileUnchecked("ZEMvsZDC"));
 133   }
 134   AliEventplane *ep = dynamic_cast<AliEventplane*>(objsout->FindObject("Eventplane"));
 135   if (ep) {
 136     ep->Reset();
 137     AliEventplane *epin = esdin->GetEventplane();
 138     if (!fDoSaveBytes) {
 139       *ep = *epin;
 140     } else {
 141       if (epin->GetQVector()) {
 142         ep->SetQVector(new TVector2(*epin->GetQVector()));
 143         ep->SetEventplaneQ(epin->GetEventplane("Q"));
 144         ep->SetQsub(new TVector2(*epin->GetQsub1()),new TVector2(*epin->GetQsub2()));
 145         ep->SetQsubRes(epin->GetQsubRes());
 146       }
 147     }
 148   }
 149   if (fDoZDC) {
 150     AliESDZDC *zdc = dynamic_cast<AliESDZDC*>(objsout->FindObject("AliESDZDC"));
 151     if (zdc) {
 152       am->LoadBranch("AliESDZDC.");
 153       *zdc = *esdin->GetESDZDC();
 154     }
 155   }
 156   if (fDoV0) {
 157     AliESDVZERO *v0 = dynamic_cast<AliESDVZERO*>(objsout->FindObject("AliESDVZERO"));
 158     if (v0) {
 159       am->LoadBranch("AliESDVZERO.");
 160       *v0 = *esdin->GetVZEROData();
 161     }
 162   }
 163   if (fDoT0) {
 164     AliESDTZERO *t0 = dynamic_cast<AliESDTZERO*>(objsout->FindObject("AliESDTZERO"));
 165     if (t0) {
 166       am->LoadBranch("AliESDTZERO.");
 167       *t0 = *esdin->GetESDTZERO();
 168     }
 169   }
 170   if (fDoTPCv) {
 171     AliESDVertex *tpcv = dynamic_cast<AliESDVertex*>(objsout->FindObject("TPCVertex"));
 172     if (tpcv) {
 173       am->LoadBranch("TPCVertex.");
 174       *tpcv = *esdin->GetPrimaryVertexTPC();
 175     }
 176   }
 177   if (fDoSPDv) {
 178     AliESDVertex *spdv = dynamic_cast<AliESDVertex*>(objsout->FindObject("SPDVertex"));
 179     if (spdv) {
 180       am->LoadBranch("SPDVertex.");
 181       *spdv = *esdin->GetPrimaryVertexSPD();
 182     }
 183   }
 184   if (fDoPriv) {
 185     AliESDVertex *priv = dynamic_cast<AliESDVertex*>(objsout->FindObject("PrimaryVertex"));
 186     if (priv) {
 187       am->LoadBranch("PrimaryVertex.");
 188       *priv = *esdin->GetPrimaryVertexTracks();
 189     }
 190   }
 191   if (fDoEmCs) {
 192     AliESDCaloCells *ecells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("EMCALCells"));
 193     if (ecells) {
 194       am->LoadBranch("EMCALCells.");
 195       *ecells = *esdin->GetEMCALCells();
 196     }
 197   }
 198   if (fDoPCs) {
 199     AliESDCaloCells *pcells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("PHOSCells"));
 200     if (pcells) {
 201       am->LoadBranch("PHOSCells.");
 202       *pcells = *esdin->GetPHOSCells();
 203     }
 204   }
 205   if (fDoEmT) {
 206     AliESDCaloTrigger *etrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("EMCALTrigger"));
 207     if (etrig) {
 208       am->LoadBranch("EMCALTrigger.");
 209       *etrig = *esdin->GetCaloTrigger("EMCAL");
 210     }
 211   }
 212   if (fDoPT) {
 213     AliESDCaloTrigger *ptrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("PHOSTrigger"));
 214     if (ptrig) {
 215       am->LoadBranch("PHOSTrigger.");
 216       *ptrig = *esdin->GetCaloTrigger("PHOS");
 217     }
 218   }
 219   if (fDoFmd) {
 220     AliESDFMD *fmd = dynamic_cast<AliESDFMD*>(objsout->FindObject("AliESDFMD"));
 221     if (fmd) {
 222       am->LoadBranch("AliESDFMD.");
 223       if (!fDoSaveBytes) {
 224         *fmd = *esdin->GetFMDData();
 225       }
 226     }
 227   }
 228   if (fDoMult) {
 229     AliMultiplicity *mult = dynamic_cast<AliMultiplicity*>(objsout->FindObject("AliMultiplicity"));
 230     if (mult) {
 231       am->LoadBranch("AliMultiplicity.");
 232       if (!fDoSaveBytes) {
 233         *mult = *esdin->GetMultiplicity();
 234       } else {
 235         const AliMultiplicity *multin = esdin->GetMultiplicity();;
 236         mult->SetFiredChips(0, multin->GetNumberOfFiredChips(0));
 237         mult->SetFiredChips(1, multin->GetNumberOfFiredChips(1));
 238         for (Int_t i=0; i<6; ++i)
 239           mult->SetITSClusters(i,multin->GetNumberOfITSClusters(i));
 240       }
 241     }
 242   }
 243   if (fDoTof) {
 244     AliTOFHeader *tofh = dynamic_cast<AliTOFHeader*>(objsout->FindObject("AliTOFHeader"));
 245     if (tofh) {
 246       am->LoadBranch("AliTOFHeader.");
 247       *tofh = *esdin->GetTOFHeader();
 248     }
 249   }
 250   if (fDoPileup) {
 251     TClonesArray *spup = dynamic_cast<TClonesArray*>(objsout->FindObject("SPDPileupVertices"));
 252     if (spup) {
 253       am->LoadBranch("SPDPileupVertices");
 254       Int_t N = esdin->GetNumberOfPileupVerticesSPD();
 255       for (Int_t i=0; i<N; ++i) {
 256         const AliESDVertex *vtx = esdin->GetPileupVertexSPD(i);
 257         if (vtx)
 258           fEvent->AddPileupVertexSPD(vtx);
 259       }
 260     }
 261     TClonesArray *tpup = dynamic_cast<TClonesArray*>(objsout->FindObject("TrkPileupVertices"));
 262     if (tpup) {
 263       am->LoadBranch("TrkPileupVertices");
 264       Int_t N = esdin->GetNumberOfPileupVerticesTracks();
 265       for (Int_t i=0; i<N; ++i) {
 266         const AliESDVertex *vtx = esdin->GetPileupVertexTracks(i);
 267         if (vtx)
 268           fEvent->AddPileupVertexTracks(vtx);
 269       }
 270     }
 271   }
 272   if (fDoClus) {
 273     TClonesArray *clus = dynamic_cast<TClonesArray*>(objsout->FindObject("CaloClusters"));
 274     if (clus) {
 275       am->LoadBranch("CaloClusters");
 276       Int_t N = esdin->GetNumberOfCaloClusters();
 277       for (Int_t i=0; i<N; ++i) {
 278         AliESDCaloCluster *c = esdin->GetCaloCluster(i);
 279         if (!c)
 280           continue;
 281         if (fPhosClusOnly && c->IsEMCAL())
 282           continue;
 283         if (fEmcalClusOnly && c->IsPHOS())
 284           continue;
 285         fEvent->AddCaloCluster(c);
 286       }
 287     }
 288   }
 289   TObjArray *namearr = fEmcNames.Tokenize(";");
 290   if (namearr) {
 291     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
 292       TString cname(namearr->At(i)->GetName());
 293       if (cname.Length()<=0)
 294         continue;
 295       TClonesArray *arrin  = dynamic_cast<TClonesArray*>(objsin->FindObject(cname));
 296       if (!arrin) {
 297         AliFatal(Form("Can not find input clusters with name %s", cname.Data()));
 298         return;
 299       }
 300       TClonesArray *arrout = dynamic_cast<TClonesArray*>(objsout->FindObject(cname));
 301       if (!arrout) {
 302         AliFatal(Form("Can not find output clusters with name %s", cname.Data()));
 303         return;
 304       }
 305       arrout->Delete();
 306       Double_t emin=0.1;
 307       if (cname.Contains("FEE"))
 308         emin = 1;
 309       const Int_t N = arrin->GetEntries();
 310       for (Int_t iC=0, nC=0; iC<N; ++iC) {
 311         AliESDCaloCluster *c = dynamic_cast<AliESDCaloCluster*>(arrin->At(iC));
 312         if (!c)
 313           continue;
 314         if (c->E()<emin)
 315           continue;
 316         AliESDCaloCluster *newCluster = new ((*arrout)[nC]) AliESDCaloCluster(*c);
 317         newCluster->SetID(nC);
 318         ++nC;
 319       }
 320     }
 321     delete namearr;
 322   }
 323   if (fDoTracks) {
 324     if (fTracks == "Tracks")
 325       am->LoadBranch("Tracks");
 326     TClonesArray *tracksin = dynamic_cast<TClonesArray*>(objsin->FindObject(fTracks));
 327     if (!tracksin) {
 328       AliFatal(Form("Can not find tracks with name %s", fTracks.Data()));
 329       return;
 330     }
 331     TString trkoutname("Tracks");
 332     if (fDoPicoTracks)
 333       trkoutname = "PicoTracks";
 334     TClonesArray *tracksout = 0;
 335     tracksout = dynamic_cast<TClonesArray*>(objsout->FindObject(trkoutname));
 336     if (!tracksout) {
 337       AliFatal(Form("Can not find tracks with name %s", trkoutname.Data()));
 338       return;
 339     }
 340     tracksout->Delete();
 341     TString classname(tracksin->GetClass()->GetName());
 342     if (classname == "AliPicoTrack") {
 343       if (!fDoPicoTracks) {
 344         AliFatal(Form("Need to enable fDoPicoTracks when reading AliPicoTracks"));
 345         return;
 346       }
 347       const Int_t Ntracks = tracksin->GetEntries();
 348       for (Int_t iTracks = 0, nacc = 0; iTracks < Ntracks; ++iTracks) {
 349         AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(tracksin->At(iTracks));
 350         if (!track)
 351           continue;
 352         new ((*tracksout)[nacc]) AliPicoTrack(*track);
 353         ++nacc;
 354       }
 355     } else {
 356       const Int_t Ntracks = tracksin->GetEntries();
 357       for (Int_t iTracks = 0, nacc = 0; iTracks < Ntracks; ++iTracks) {
 358         AliESDtrack *track = dynamic_cast<AliESDtrack*>(tracksin->At(iTracks));
 359         if (!track)
 360           continue;
 361         if (fCuts) {
 362           if (!fCuts->IsSelected(track))
 363             continue;
 364         }
 365         if (fDoPicoTracks) {
 366           AliEsdTrackExt newtrack(*track);
 367           Double_t etaemc = 0;
 368           Double_t phiemc = 0;
 369           if (newtrack.IsEMCAL()) {
 370             etaemc = newtrack.GetEmcEta();
 371             phiemc = newtrack.GetEmcPhi();
 372           }
 373           new ((*tracksout)[nacc]) AliPicoTrack(newtrack.Pt(), newtrack.Eta(), newtrack.Phi(),
 374                                                 newtrack.Charge(), newtrack.GetLabel(),
 375                                                 etaemc, phiemc, newtrack.IsEMCAL());
 376           ++nacc;
 377         } else {
 378           AliEsdTrackExt *newtrack = new ((*tracksout)[nacc]) AliEsdTrackExt(*track);
 379           if (fDoMiniTracks) {
 380             newtrack->MakeMiniTrack(0,fRemoveCP);
 381             newtrack->ResetCovariance(fResetCov);
 382           } else {
 383             newtrack->DeleteParams();
 384           }
 385           newtrack->SetID(nacc);
 386           ++nacc;
 387         }
 388       }
 389     }
 390   }
 391   if (fDoMuonTracks) {
 392     TClonesArray *muons = dynamic_cast<TClonesArray*>(objsout->FindObject("MuonTracks"));
 393     if (muons) {
 394       am->LoadBranch("MuonTracks");
 395       Int_t N = esdin->GetNumberOfMuonTracks();
 396       for (Int_t i=0; i<N; ++i) {
 397         const AliESDMuonTrack *mt = esdin->GetMuonTrack(i);
 398         new ((*muons)[i]) AliESDMuonTrack(*mt);
 399       }
 400     }
 401   }
 402   fTree->Fill();
 403 }
 404
 405 //_________________________________________________________________________________________________
 406 void AliEsdSkimTask::UserCreateOutputObjects()
 407 {
 408   // Create output objects.
 409
 410   TFile *file = OpenFile(1);
 411   fTree = new TTree("esdTree", "Tree with skimmed ESD objects");
 412   file->SetCompressionLevel(5);
 413   fTree->SetDirectory(file);
 414   fTree->SetAutoFlush(-10*1024*1024);
 415
 416   fEvent = new AliESDEvent;
 417   fEvent->CreateStdContent();
 418   if (fDoCent) {
 419     AliCentrality *cent = new AliCentrality;
 420     fEvent->AddObject(cent);
 421   }
 422   if (fDoRP) {
 423     AliEventplane *ep = new AliEventplane;
 424     fEvent->AddObject(ep);
 425   }
 426   TObjArray *namearr = fEmcNames.Tokenize(";");
 427   if (namearr) {
 428     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
 429       TString cname(namearr->At(i)->GetName());
 430       if (cname.Length()<=0)
 431         continue;
 432       TClonesArray *arr = new TClonesArray("AliESDCaloCluster",0);
 433       arr->SetName(cname);
 434       fEvent->AddObject(arr);
 435     }
 436     delete namearr;
 437   }
 438   if (fDoTracks) {
 439     if (fDoPicoTracks) {
 440       TClonesArray *arr = new TClonesArray("AliPicoTrack",0);
 441       arr->SetName("PicoTracks");
 442       fEvent->AddObject(arr);
 443     } else {
 444       TObject *obj = fEvent->GetList()->FindObject("Tracks");
 445       fEvent->GetList()->Remove(obj);
 446       TClonesArray *arr = new TClonesArray("AliEsdTrackExt",0);
 447       arr->SetName("Tracks");
 448       fEvent->AddObject(arr);
 449     }
 450   }
 451   fEvent->GetStdContent();
 452   fEvent->WriteToTree(fTree);
 453   fTree->GetUserInfo()->Add(fEvent);
 454   PostData(1,fTree);
 455 }