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