PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAflow.cxx

   1 #include "TMath.h"
   2 #include "TH1D.h"
   3 #include "TH2D.h"
   4 #include "TSeqCollection.h"
   5 #include "TObjArray.h"
   6 #include "TObjArray.h"
   7 #include "TChain.h"
   8 #include "TMCProcess.h"
   9 #include "TLorentzVector.h"
  10 #include "TDirectory.h"
  11 #include "TROOT.h"
  12
  13 #include "AliLog.h"
  14 #include "AliVParticle.h"
  15 #include "AliMCParticle.h"
  16 #include "AliStack.h"
  17 #include "AliVEvent.h"
  18 #include "AliESDEvent.h"
  19 #include "AliMCEvent.h"
  20 #include "AliESDtrack.h"
  21 #include "AliFlowTrackCuts.h"
  22 #include "AliFlowEventCuts.h"
  23 #include "AliMultiplicity.h"
  24 #include "AliESDtrackCuts.h"
  25 #include "AliVertex.h"
  26
  27 #include "AliAnalysisTask.h"
  28 #include "AliAnalysisManager.h"
  29
  30 #include "AliAnalysisTaskQAflow.h"
  31
  32 ClassImp(AliAnalysisTaskQAflow)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow()
  36   : AliAnalysisTaskSE(),
  37     fOutput(NULL),
  38     fEventCuts(NULL),
  39     fTrackCuts(NULL)
  40 {
  41   // Default constructor
  42 }
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow(const char* name)
  46   : AliAnalysisTaskSE(name),
  47     fOutput(NULL),
  48     fEventCuts(NULL),
  49     fTrackCuts(NULL)
  50 {
  51   // Constructor
  52   DefineOutput(1, TObjArray::Class());
  53 }
  54
  55 //________________________________________________________________________
  56 void AliAnalysisTaskQAflow::UserCreateOutputObjects()
  57 {
  58   // Called once at the beginning
  59   fOutput=new TObjArray();
  60   //TDirectory* before = gDirectory->mkdir("before cuts","before cuts");
  61   //TDirectory* after = gDirectory->mkdir("after cuts","after cuts");
  62   TObjArray* before = new TObjArray();
  63   TObjArray* after = new TObjArray();
  64   fOutput->Add(before);
  65   fOutput->Add(after);
  66
  67   //define histograms
  68
  69   TH1* hist;
  70   hist = new TH1I("tracklet_multiplicity","tracklet multiplicity",1000,0,20000);
  71   before->Add(hist); after->Add(hist->Clone()); //0
  72   hist = new TH1I("track_multiplicity", "standard TPC track multiplicity",1000,0,10000);
  73   before->Add(hist); after->Add(hist->Clone()); //1
  74   hist = new TH1I("refmult","refmult",1000, 0,10000);
  75   before->Add(hist); after->Add(hist->Clone()); //2
  76   hist = new TH1I("SPD_clusters","SPD clusters",1000,0,100000);
  77   before->Add(hist); after->Add(hist->Clone()); //3
  78   hist = new TH1D("primary_vertexZ","primary vertex z",100,-20,20);
  79   before->Add(hist); after->Add(hist->Clone()); //4
  80   hist = new TH1I("ITS_clusters_on_track", "ITS clusters on track", 8, 0, 8);
  81   before->Add(hist); after->Add(hist->Clone()); //5
  82   hist = new TH1I("TPC_clusters_on_track", "TPC clusters on track", 159, 1, 160);
  83   before->Add(hist); after->Add(hist->Clone()); //6
  84   hist = new TH1D("TPC_chi2_per_cluster","TPC #chi^{2}/cluster",100,0.0,5.0);
  85   before->Add(hist); after->Add(hist->Clone()); //7
  86   hist = new TH1D("DCA_xy","DCA xy", 1000, -5.0, 5.0 );
  87   before->Add(hist); after->Add(hist->Clone()); //8
  88   hist = new TH1D("DCA_z","DCA z", 1000, -5.0, 5.0 );
  89   before->Add(hist); after->Add(hist->Clone()); //9
  90   hist = new TH1D("phi_tracklets","#phi tracklets", 1000, 0.0, TMath::TwoPi() );
  91   before->Add(hist); after->Add(hist->Clone()); //10
  92   hist = new TH1D("phi_tracks","#phi tracks", 1000, 0.0, TMath::TwoPi() );
  93   before->Add(hist); after->Add(hist->Clone()); //11
  94   hist = new TH1D("eta_tracklets","#eta tracklets", 1000, -2.0, 2.0 );
  95   before->Add(hist); after->Add(hist->Clone()); //12
  96   hist = new TH1D("eta_tracks","#eta tracks", 1000, -2.0, 2.0 );
  97   before->Add(hist); after->Add(hist->Clone()); //13
  98   hist = new TH1D("TPC_vertex_z", "TPC vertex z", 100,-20.0,20.0);
  99   before->Add(hist); after->Add(hist->Clone()); //14
 100
 101   //post data here as it doesn't change anyway (the pointer to list anyway)
 102
 103   //restore dir add status
 104   PostData(1, fOutput);
 105 }
 106
 107 //________________________________________________________________________
 108 void AliAnalysisTaskQAflow::UserExec(Option_t *)
 109 {
 110
 111   //get teh input data
 112   AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
 113   if (!event)
 114   {
 115     AliFatal("no ESD event");
 116     return;
 117   }
 118
 119   //TObjArray* before = ((TDirectory*)fOutput->At(0))->GetList();
 120   //TObjArray* after = ((TDirectory*)fOutput->At(1))->GetList();
 121   TObjArray* before = (TObjArray*)fOutput->At(0);
 122   TObjArray* after = (TObjArray*)fOutput->At(1);
 123   TH1* htrackletmultB = dynamic_cast<TH1*>(before->At(0));
 124   TH1* htrackletmultA = dynamic_cast<TH1*>(after->At(0));
 125   TH1* htrackmultB = dynamic_cast<TH1*>(before->At(1));
 126   TH1* htrackmultA = dynamic_cast<TH1*>(after->At(1));
 127   TH1* hrefmultB = dynamic_cast<TH1*>(before->At(2));
 128   TH1* hrefmultA = dynamic_cast<TH1*>(after->At(2));
 129   TH1* hspdclustersB = dynamic_cast<TH1*>(before->At(3));
 130   TH1* hspdclustersA = dynamic_cast<TH1*>(after->At(3));
 131   TH1* hprimvtxzB = dynamic_cast<TH1*>(before->At(4));
 132   TH1* hprimvtxzA = dynamic_cast<TH1*>(after->At(4));
 133
 134   TH1* hITSclsB = dynamic_cast<TH1*>(before->At(5));
 135   TH1* hITSclsA = dynamic_cast<TH1*>(after->At(5));
 136   TH1* hTPCclsB = dynamic_cast<TH1*>(before->At(6));
 137   TH1* hTPCclsA = dynamic_cast<TH1*>(after->At(6));
 138   TH1* hTPCchi2B = dynamic_cast<TH1*>(before->At(7));
 139   TH1* hTPCchi2A = dynamic_cast<TH1*>(after->At(7));
 140   TH1* hdcaxyB = dynamic_cast<TH1*>(before->At(8));
 141   TH1* hdcaxyA = dynamic_cast<TH1*>(after->At(8));
 142   TH1* hdcazB = dynamic_cast<TH1*>(before->At(9));
 143   TH1* hdcazA = dynamic_cast<TH1*>(after->At(9));
 144   TH1* hphitrackletsB = dynamic_cast<TH1*>(before->At(10));
 145   TH1* hphitrackletsA = dynamic_cast<TH1*>(after->At(10));
 146   TH1* hphitracksB = dynamic_cast<TH1*>(before->At(11));
 147   TH1* hphitracksA = dynamic_cast<TH1*>(after->At(11));
 148   TH1* hetatrackletsB = dynamic_cast<TH1*>(before->At(12));
 149   TH1* hetatrackletsA = dynamic_cast<TH1*>(after->At(12));
 150   TH1* hetatracksB = dynamic_cast<TH1*>(before->At(13));
 151   TH1* hetatracksA = dynamic_cast<TH1*>(after->At(13));
 152   TH1* hprimvtxzTPCB = dynamic_cast<TH1*>(before->At(14));
 153   TH1* hprimvtxzTPCA = dynamic_cast<TH1*>(after->At(14));
 154
 155   AliMultiplicity* tracklets = const_cast<AliMultiplicity*>(event->GetMultiplicity());
 156   Int_t ntracklets=0;
 157   Int_t nspdclusters=0;
 158   if (tracklets)
 159   {
 160     ntracklets = tracklets->GetNumberOfTracklets();
 161     nspdclusters = tracklets->GetNumberOfITSClusters(0,1);
 162     for (Int_t i=0; i<tracklets->GetNumberOfTracklets(); i++)
 163     {
 164       Bool_t pass=fTrackCuts->IsSelected(tracklets,i);
 165       Float_t phi=tracklets->GetPhi(i);
 166       Float_t eta=tracklets->GetEta(i);
 167       hphitrackletsB->Fill(phi); if (pass) hphitrackletsA->Fill(phi);
 168       hetatrackletsB->Fill(eta); if (pass) hetatrackletsA->Fill(eta);
 169     }
 170
 171   }
 172   Int_t trackmult = AliESDtrackCuts::GetReferenceMultiplicity(event,kTRUE);
 173   Int_t refmult = fEventCuts->RefMult(event);
 174   Bool_t passevent = fEventCuts->IsSelected(event);
 175   htrackletmultB->Fill(ntracklets); if (passevent) htrackletmultA->Fill(ntracklets);
 176   htrackmultB->Fill(trackmult); if (passevent) htrackmultA->Fill( trackmult);
 177   hrefmultB->Fill(refmult); if (passevent) hrefmultA->Fill( refmult);
 178   hspdclustersB->Fill(nspdclusters); if (passevent) hspdclustersA->Fill( nspdclusters);
 179   const AliVertex* vertex = event->GetPrimaryVertex();
 180   Float_t vtxz=0.;
 181   if (vertex)
 182   {
 183     vtxz = vertex->GetZ();
 184     hprimvtxzB->Fill(vtxz); if (passevent) hprimvtxzA->Fill(vtxz);
 185   }
 186   const AliVertex* vertextpc = event->GetPrimaryVertexTPC();
 187   if (vertextpc)
 188   {
 189     vtxz = vertextpc->GetZ();
 190     hprimvtxzTPCB->Fill(vtxz); if (passevent) hprimvtxzTPCA->Fill(vtxz);
 191   }
 192   fTrackCuts->SetEvent(event);
 193   for (Int_t i=0; i<fTrackCuts->GetNumberOfInputObjects(); i++)
 194   {
 195     TObject* obj = fTrackCuts->GetInputObject(i);
 196     Bool_t pass = fTrackCuts->IsSelected(obj,i);
 197     Float_t dcaxy=0.0;
 198     Float_t dcaz=0.0;
 199     Float_t tpcchi2=0.0;
 200     Float_t tpcchi2percls=0.0;
 201     Int_t ntpccls=0;
 202     Int_t nitscls=0;
 203     Float_t eta=0.0;
 204     Float_t phi=0.0;
 205     AliESDtrack* track = dynamic_cast<AliESDtrack*>(fTrackCuts->GetTrack());
 206     if (track)
 207     {
 208       track->GetImpactParameters(dcaxy,dcaz);
 209       tpcchi2=track->GetTPCchi2();
 210       ntpccls=track->GetTPCNcls();
 211       eta=track->Eta();
 212       phi=track->Phi();
 213       tpcchi2percls= (ntpccls==0)?0.0:tpcchi2/ntpccls;
 214       nitscls = track->GetNcls(0);
 215       hITSclsB->Fill(nitscls); if (pass) hITSclsA->Fill( nitscls);
 216       hTPCclsB->Fill(ntpccls); if (pass) hTPCclsA->Fill( ntpccls);
 217       hTPCchi2B->Fill(tpcchi2percls); if (pass) hTPCchi2A->Fill( tpcchi2percls);
 218       hdcaxyB->Fill(dcaxy); if (pass) hdcaxyA->Fill( dcaxy);
 219       hdcazB->Fill(dcaz); if (pass) hdcazA->Fill(dcaz);
 220       hetatracksB->Fill(eta); if (pass) hetatracksA->Fill(eta);
 221       hphitracksB->Fill(phi); if (pass) hphitracksA->Fill(phi);
 222     }
 223   }
 224 }
 225
 226 //________________________________________________________________________
 227 void AliAnalysisTaskQAflow::Terminate(Option_t *)
 228 {
 229
 230 }
 231
 232 //________________________________________________________________________
 233 AliAnalysisTaskQAflow::~AliAnalysisTaskQAflow()
 234 {
 235   //dtor
 236   delete fTrackCuts;
 237   delete fEventCuts;
 238 }