PWG4/UserTasks/EmcalTasks/AliEsdSkimTask.cxx

   1 // $Id$
   2 //
   3 // Task to skim ESD files.
   4 //
   5 //
   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 "AliEventplane.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDFMD.h"
  16 #include "AliESDtrackCuts.h"
  17 #include "AliEsdTrackExt.h"
  18 #include "AliMultiplicity.h"
  19
  20 //_________________________________________________________________________________________________
  21 AliEsdSkimTask::AliEsdSkimTask(const char *opt) :
  22   AliAnalysisTaskSE(opt), fEvent(0), fTree(0), fCuts(0),
  23   fDoZDC(1), fDoV0(1), fDoT0(1), fDoTPCv(1), fDoSPDv(1), fDoPriv(1),
  24   fDoEmCs(1), fDoPCs(1), fDoEmT(1), fDoPT(1), fDoTracks(1), fDoFmd(1),
  25   fDoMult(1), fDoTof(1), fDoPileup(1), fDoClus(1), fEmcNames(""),
  26   fDoMiniTracks(0), fTracks("Tracks"), fPhosClusOnly(0), fDoSaveBytes(1),
  27   fDoCent(1), fDoRP(1)
  28 {
  29   // Constructor.
  30
  31   if (!opt)
  32     return;
  33
  34   fBranchNames = "ESD:AliESDHeader.,AliESDRun.";
  35
  36   DefineOutput(1, TTree::Class());
  37 }
  38
  39 //_________________________________________________________________________________________________
  40 void AliEsdSkimTask::UserExec(Option_t */*opt*/)
  41 {
  42   // Process event.
  43
  44   AliESDEvent *esdin = dynamic_cast<AliESDEvent*>(InputEvent());
  45   if (!esdin)
  46     return;
  47
  48   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  49
  50   fEvent->Reset();
  51
  52   TList* objsin  = esdin->GetList();
  53   TList* objsout = fEvent->GetList();
  54
  55   AliESDHeader *header = dynamic_cast<AliESDHeader*>(objsout->FindObject("AliESDHeader"));
  56   if (header) {
  57     am->LoadBranch("AliESDHeader.");
  58     *header = *esdin->GetHeader();
  59   }
  60   AliESDRun *run = dynamic_cast<AliESDRun*>(objsout->FindObject("AliESDRun"));
  61   if (run) {
  62     am->LoadBranch("AliESDRun.");
  63     *run = *esdin->GetESDRun();
  64   }
  65   AliCentrality *cent = dynamic_cast<AliCentrality*>(objsout->FindObject("Cent"));
  66   if (cent) {
  67     cent->Reset();
  68     AliCentrality *centin = esdin->GetCentrality();
  69     cent->SetQuality(centin->GetQuality());
  70     cent->SetCentralityV0M(centin->GetCentralityPercentileUnchecked("V0M"));
  71     cent->SetCentralityFMD(centin->GetCentralityPercentileUnchecked("FMD"));
  72     cent->SetCentralityTRK(centin->GetCentralityPercentileUnchecked("TRK"));
  73     cent->SetCentralityTKL(centin->GetCentralityPercentileUnchecked("TKL"));
  74     cent->SetCentralityCL0(centin->GetCentralityPercentileUnchecked("CL0"));
  75     cent->SetCentralityCL1(centin->GetCentralityPercentileUnchecked("CL1"));
  76     cent->SetCentralityV0MvsFMD(centin->GetCentralityPercentileUnchecked("V0MvsFMD"));
  77     cent->SetCentralityTKLvsV0M(centin->GetCentralityPercentileUnchecked("TKLvsV0M"));
  78     cent->SetCentralityZEMvsZDC(centin->GetCentralityPercentileUnchecked("ZEMvsZDC"));
  79   }
  80   AliEventplane *ep = dynamic_cast<AliEventplane*>(objsout->FindObject("EP"));
  81   if (ep) {
  82     ep->Reset();
  83     AliEventplane *epin = esdin->GetEventplane();
  84     if (!fDoSaveBytes) {
  85       *ep = *epin;
  86     } else {
  87       if (epin->GetQVector()) {
  88         ep->SetQVector(new TVector2(*epin->GetQVector()));
  89         ep->SetEventplaneQ(epin->GetEventplane("Q"));
  90         ep->SetQsub(new TVector2(*epin->GetQsub1()),new TVector2(*epin->GetQsub2()));
  91         ep->SetQsubRes(epin->GetQsubRes());
  92       }
  93     }
  94   }
  95   AliESDZDC *zdc = dynamic_cast<AliESDZDC*>(objsout->FindObject("AliESDZDC"));
  96   if (zdc) {
  97     am->LoadBranch("AliESDZDC.");
  98     *zdc = *esdin->GetESDZDC();
  99   }
 100   AliESDVZERO *v0 = dynamic_cast<AliESDVZERO*>(objsout->FindObject("AliESDVZERO"));
 101   if (v0) {
 102     am->LoadBranch("AliESDVZERO.");
 103     *v0 = *esdin->GetVZEROData();
 104   }
 105   AliESDTZERO *t0 = dynamic_cast<AliESDTZERO*>(objsout->FindObject("AliESDTZERO"));
 106   if (t0) {
 107     am->LoadBranch("AliESDTZERO.");
 108     *t0 = *esdin->GetESDTZERO();
 109   }
 110   AliESDVertex *tpcv = dynamic_cast<AliESDVertex*>(objsout->FindObject("TPCVertex"));
 111   if (tpcv) {
 112     am->LoadBranch("TPCVertex.");
 113     *tpcv = *esdin->GetPrimaryVertexTPC();
 114   }
 115   AliESDVertex *spdv = dynamic_cast<AliESDVertex*>(objsout->FindObject("SPDVertex"));
 116   if (spdv) {
 117     am->LoadBranch("SPDVertex.");
 118     *spdv = *esdin->GetPrimaryVertexSPD();
 119   }
 120   AliESDVertex *priv = dynamic_cast<AliESDVertex*>(objsout->FindObject("PrimaryVertex"));
 121   if (priv) {
 122     am->LoadBranch("PrimaryVertex.");
 123     *priv = *esdin->GetPrimaryVertexTracks();
 124   }
 125   AliESDCaloCells *ecells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("EMCALCells"));
 126   if (ecells) {
 127     am->LoadBranch("EMCALCells.");
 128     *ecells = *esdin->GetEMCALCells();
 129   }
 130   AliESDCaloCells *pcells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("PHOSCells"));
 131   if (pcells) {
 132     am->LoadBranch("PHOSCells.");
 133     *pcells = *esdin->GetPHOSCells();
 134   }
 135   AliESDCaloTrigger *etrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("EMCALTrigger"));
 136   if (etrig) {
 137     am->LoadBranch("EMCALTrigger.");
 138     *etrig = *esdin->GetCaloTrigger("EMCAL");
 139   }
 140   AliESDCaloTrigger *ptrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("PHOSTrigger"));
 141   if (ptrig) {
 142     am->LoadBranch("PHOSTrigger.");
 143     *ptrig = *esdin->GetCaloTrigger("PHOS");
 144   }
 145   AliESDFMD *fmd = dynamic_cast<AliESDFMD*>(objsout->FindObject("AliESDFMD"));
 146   if (fmd) {
 147     am->LoadBranch("AliESDFMD.");
 148     if (!fDoSaveBytes) {
 149       *fmd = *esdin->GetFMDData();
 150     }
 151   }
 152   AliMultiplicity *mult = dynamic_cast<AliMultiplicity*>(objsout->FindObject("AliMultiplicity"));
 153   if (mult) {
 154     am->LoadBranch("AliMultiplicity.");
 155     if (!fDoSaveBytes) {
 156       *mult = *esdin->GetMultiplicity();
 157     } else {
 158       const AliMultiplicity *multin = esdin->GetMultiplicity();;
 159       mult->SetFiredChips(0, multin->GetNumberOfFiredChips(0));
 160       mult->SetFiredChips(1, multin->GetNumberOfFiredChips(1));
 161       for (Int_t i=0; i<6; ++i)
 162         mult->SetITSClusters(i,mult->GetNumberOfITSClusters(i));
 163     }
 164   }
 165   AliTOFHeader *tofh = dynamic_cast<AliTOFHeader*>(objsout->FindObject("AliTOFHeader"));
 166   if (tofh) {
 167     am->LoadBranch("AliTOFHeader.");
 168     *tofh = *esdin->GetTOFHeader();
 169   }
 170   TClonesArray *spup = dynamic_cast<TClonesArray*>(objsout->FindObject("SPDPileupVertices"));
 171   if (spup) {
 172     am->LoadBranch("SPDPileupVertices");
 173     Int_t N = esdin->GetNumberOfPileupVerticesSPD();
 174     for (Int_t i=0; i<N; ++i) {
 175       const AliESDVertex *vtx = esdin->GetPileupVertexSPD(i);
 176       if (vtx)
 177         fEvent->AddPileupVertexSPD(vtx);
 178     }
 179   }
 180   TClonesArray *tpup = dynamic_cast<TClonesArray*>(objsout->FindObject("TrkPileupVertices"));
 181   if (tpup) {
 182     am->LoadBranch("TrkPileupVertices");
 183     Int_t N = esdin->GetNumberOfPileupVerticesTracks();
 184     for (Int_t i=0; i<N; ++i) {
 185       const AliESDVertex *vtx = esdin->GetPileupVertexTracks(i);
 186       if (vtx)
 187         fEvent->AddPileupVertexTracks(vtx);
 188     }
 189   }
 190   TClonesArray *clus = dynamic_cast<TClonesArray*>(objsout->FindObject("CaloClusters"));
 191   if (clus) {
 192     am->LoadBranch("CaloClusters");
 193     Int_t N = esdin->GetNumberOfCaloClusters();
 194     for (Int_t i=0; i<N; ++i) {
 195       AliESDCaloCluster *c = esdin->GetCaloCluster(i);
 196       if (fPhosClusOnly && c->IsEMCAL())
 197         continue;
 198       if (c)
 199         fEvent->AddCaloCluster(c);
 200     }
 201   }
 202   TObjArray *namearr = fEmcNames.Tokenize(";");
 203   if (namearr) {
 204     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
 205       TString cname(namearr->At(i)->GetName());
 206       if (cname.Length()<=0)
 207         continue;
 208       TClonesArray *arrin  = dynamic_cast<TClonesArray*>(objsin->FindObject(cname));
 209       if (!arrin) {
 210         AliFatal(Form("Can not find input clusters with name %s", cname.Data()));
 211         return;
 212       }
 213       TClonesArray *arrout = dynamic_cast<TClonesArray*>(objsout->FindObject(cname));
 214       if (!arrout) {
 215         AliFatal(Form("Can not find output clusters with name %s", cname.Data()));
 216         return;
 217       }
 218       arrout->Delete();
 219       const Int_t N = arrin->GetEntries();
 220       for (Int_t iC=0, nC=0; iC<N; ++iC) {
 221         AliESDCaloCluster *c = dynamic_cast<AliESDCaloCluster*>(arrin->At(iC));
 222         if (!c)
 223           continue;
 224         new ((*arrout)[nC++]) AliESDCaloCluster(*c);
 225       }
 226     }
 227     delete namearr;
 228   }
 229   if (fDoTracks) {
 230     am->LoadBranch("Tracks");
 231     TClonesArray *tracksin = dynamic_cast<TClonesArray*>(objsin->FindObject(fTracks));
 232     if (!tracksin) {
 233       AliFatal(Form("Can not find tracks with name %s", fTracks.Data()));
 234       return;
 235     }
 236     TClonesArray *tracksout = dynamic_cast<TClonesArray*>(objsout->FindObject("Tracks"));
 237     if (!tracksout) {
 238       AliFatal(Form("Can not find tracks with name %s", "Tracks"));
 239       return;
 240     }
 241     const Int_t Ntracks = tracksin->GetEntries();
 242     Int_t nacc = 0;
 243     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 244       AliESDtrack *track = dynamic_cast<AliESDtrack*>(tracksin->At(iTracks));
 245       if (!track)
 246         continue;
 247       if (fCuts) {
 248         if (!fCuts->IsSelected(track))
 249           continue;
 250       }
 251       AliEsdTrackExt *newtrack = new ((*tracksout)[nacc]) AliEsdTrackExt(*track);
 252       if (fDoMiniTracks) {
 253         newtrack->MakeMiniTrack();
 254       } else {
 255         newtrack->DeleteParams();
 256       }
 257       newtrack->SetID(nacc);
 258       fEvent->AddTrack(track);
 259       ++nacc;
 260     }
 261     if (fCuts)
 262       AliInfo(Form("Selected %d out of %d \n", nacc, Ntracks));
 263   }
 264   fTree->Fill();
 265 }
 266
 267 //_________________________________________________________________________________________________
 268 void AliEsdSkimTask::UserCreateOutputObjects()
 269 {
 270   // Create output objects.
 271
 272   fTree = new TTree("esdTree", "Tree with skimmed ESD objects");
 273   fEvent = new AliESDEvent;
 274   fEvent->AddObject(new AliESDHeader());
 275   fEvent->AddObject(new AliESDRun());
 276   if (fDoCent) {
 277     AliCentrality *cent = new AliCentrality;
 278     cent->SetName("Cent");
 279     fEvent->AddObject(cent);
 280   }
 281   if (fDoRP) {
 282     AliEventplane *ep = new AliEventplane;
 283     ep->SetName("EP");
 284     fEvent->AddObject(ep);
 285   }
 286   if (fDoZDC)
 287     fEvent->AddObject(new AliESDZDC());
 288   if (fDoV0)
 289     fEvent->AddObject(new AliESDVZERO());
 290   if (fDoT0)
 291     fEvent->AddObject(new AliESDTZERO());
 292   if (fDoTPCv) {
 293     AliESDVertex *tpcv = new AliESDVertex();
 294     tpcv->SetName("TPCVertex");
 295     fEvent->AddObject(tpcv);
 296   }
 297   if (fDoSPDv) {
 298     AliESDVertex *spdv = new AliESDVertex();
 299     spdv->SetName("SPDVertex");
 300     fEvent->AddObject(spdv);
 301   }
 302   if (fDoPriv) {
 303     AliESDVertex *priv = new AliESDVertex();
 304     priv->SetName("PrimaryVertex");
 305     fEvent->AddObject(priv);
 306   }
 307   if (fDoEmCs) {
 308     fEvent->AddObject(new AliESDCaloCells("EMCALCells","EMCALCells"));
 309   }
 310   if (fDoPCs) {
 311     fEvent->AddObject(new AliESDCaloCells("PHOSCells","PHOSCells"));
 312   }
 313   if (fDoEmT) {
 314     AliESDCaloTrigger *etrig = new AliESDCaloTrigger;
 315     etrig->SetName("EMCALTrigger");
 316     fEvent->AddObject(etrig);
 317   }
 318   if (fDoPT) {
 319     AliESDCaloTrigger *ptrig = new AliESDCaloTrigger;
 320     ptrig->SetName("PHOSTrigger");
 321     fEvent->AddObject(ptrig);
 322   }
 323   if (fDoFmd) {
 324     AliESDFMD *fmd = new AliESDFMD;
 325     fEvent->AddObject(fmd);
 326   }
 327   if (fDoMult) {
 328     fEvent->AddObject(new AliMultiplicity());
 329   }
 330   if (fDoPileup) {
 331     TClonesArray *arr1 = new TClonesArray("AliESDVertex",0);
 332     arr1->SetName("SPDPileupVertices");
 333     fEvent->AddObject(arr1);
 334     TClonesArray *arr2 = new TClonesArray("AliESDVertex",0);
 335     arr2->SetName("TrkPileupVertices");
 336     fEvent->AddObject(arr2);
 337   }
 338   if (fDoTof) {
 339     fEvent->AddObject(new AliTOFHeader());
 340   }
 341   if (fDoClus) {
 342     TClonesArray *arr = new TClonesArray("AliESDCaloCluster",0);
 343     arr->SetName("CaloClusters");
 344     fEvent->AddObject(arr);
 345   }
 346   TObjArray *namearr = fEmcNames.Tokenize(";");
 347   if (namearr) {
 348     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
 349       TString cname(namearr->At(i)->GetName());
 350       if (cname.Length()<=0)
 351         continue;
 352       TClonesArray *arr = new TClonesArray("AliESDCaloCluster",0);
 353       arr->SetName(cname);
 354       fEvent->AddObject(arr);
 355     }
 356     delete namearr;
 357   }
 358   if (fDoTracks) {
 359     TClonesArray *arr = 0;
 360     if (fDoMiniTracks) {
 361       arr = new TClonesArray("AliEsdTrackExt",0);
 362       arr->SetName("Tracks");
 363     } else {
 364       arr = new TClonesArray("AliESDtrackExt",0);
 365       arr->SetName("Tracks");
 366     }
 367     fEvent->AddObject(arr);
 368   }
 369   fEvent->GetStdContent();
 370   fEvent->WriteToTree(fTree);
 371   fTree->GetUserInfo()->Add(fEvent);
 372   TFile *file = OpenFile(1);
 373   file->SetCompressionLevel(2);
 374   fTree->SetDirectory(file);
 375   fTree->SetAutoFlush(-1024*1024*1024);
 376   fTree->SetAutoSave(-1024*1024*1024);
 377   PostData(1,fTree);
 378 }