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 #include "TNtuple.h"
  13
  14 #include "AliLog.h"
  15 #include "AliVParticle.h"
  16 #include "AliMCParticle.h"
  17 #include "AliStack.h"
  18 #include "AliVEvent.h"
  19 #include "AliESDEvent.h"
  20 #include "AliMCEvent.h"
  21 #include "AliESDtrack.h"
  22 #include "AliFlowTrackCuts.h"
  23 #include "AliFlowEventCuts.h"
  24 #include "AliMultiplicity.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliVertex.h"
  27 #include "AliFlowEventSimple.h"
  28 #include "AliFlowVector.h"
  29
  30 #include "AliAnalysisTask.h"
  31 #include "AliAnalysisManager.h"
  32
  33 #include "AliAnalysisTaskQAflow.h"
  34
  35 ClassImp(AliAnalysisTaskQAflow)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow()
  39   : AliAnalysisTaskSE(),
  40     fOutput(NULL),
  41     fFillNtuple(kFALSE),
  42     fNtuple(NULL),
  43     fEventCuts(NULL),
  44     fTrackCuts(NULL)
  45 {
  46   // Default constructor
  47 }
  48
  49 //________________________________________________________________________
  50 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow(const char* name)
  51   : AliAnalysisTaskSE(name),
  52     fOutput(NULL),
  53     fFillNtuple(kFALSE),
  54     fNtuple(NULL),
  55     fEventCuts(NULL),
  56     fTrackCuts(NULL)
  57 {
  58   // Constructor
  59   DefineInput(1, AliFlowEventSimple::Class());
  60   DefineOutput(1, TObjArray::Class());
  61 }
  62
  63 //________________________________________________________________________
  64 void AliAnalysisTaskQAflow::UserCreateOutputObjects()
  65 {
  66   // Called once at the beginning
  67   fOutput=new TObjArray();
  68   fNtuple = new TNtuple("flowQAtree","flowQAtree","meanpt:qx:qy:mult:refmult:trig:tpcvtxx:tpcvtxy:tpcvtxz:ntrackletsA:ntrackletsC");
  69   //TDirectory* before = gDirectory->mkdir("before cuts","before cuts");
  70   //TDirectory* after = gDirectory->mkdir("after cuts","after cuts");
  71   TObjArray* before = new TObjArray();
  72   TObjArray* after = new TObjArray();
  73   fOutput->Add(before);
  74   fOutput->Add(after);
  75
  76   //define histograms
  77
  78   TH1* hist;
  79   hist = new TH1I("tracklet_multiplicity","tracklet multiplicity",1000,0,20000);
  80   before->Add(hist); after->Add(hist->Clone()); //0
  81   hist = new TH1I("track_multiplicity", "standard TPC track multiplicity",1000,0,10000);
  82   before->Add(hist); after->Add(hist->Clone()); //1
  83   hist = new TH1I("refmult","refmult",1000, 0,10000);
  84   before->Add(hist); after->Add(hist->Clone()); //2
  85   hist = new TH1I("SPD_clusters","SPD clusters",1000,0,100000);
  86   before->Add(hist); after->Add(hist->Clone()); //3
  87   hist = new TH1D("primary_vertexZ","primary vertex z",100,-20,20);
  88   before->Add(hist); after->Add(hist->Clone()); //4
  89   hist = new TH1I("ITS_clusters_on_track", "ITS clusters on track", 8, 0, 8);
  90   before->Add(hist); after->Add(hist->Clone()); //5
  91   hist = new TH1I("TPC_clusters_on_track", "TPC clusters on track", 159, 1, 160);
  92   before->Add(hist); after->Add(hist->Clone()); //6
  93   hist = new TH1D("TPC_chi2_per_cluster","TPC #chi^{2}/cluster",100,0.0,5.0);
  94   before->Add(hist); after->Add(hist->Clone()); //7
  95   hist = new TH1D("DCA_xy","DCA xy", 1000, -5.0, 5.0 );
  96   before->Add(hist); after->Add(hist->Clone()); //8
  97   hist = new TH1D("DCA_z","DCA z", 1000, -5.0, 5.0 );
  98   before->Add(hist); after->Add(hist->Clone()); //9
  99   hist = new TH1D("phi_tracklets","#phi tracklets", 1000, 0.0, TMath::TwoPi() );
 100   before->Add(hist); after->Add(hist->Clone()); //10
 101   hist = new TH1D("phi_tracks","#phi tracks", 1000, 0.0, TMath::TwoPi() );
 102   before->Add(hist); after->Add(hist->Clone()); //11
 103   hist = new TH1D("eta_tracklets","#eta tracklets", 1000, -2.0, 2.0 );
 104   before->Add(hist); after->Add(hist->Clone()); //12
 105   hist = new TH1D("eta_tracks","#eta tracks", 1000, -2.0, 2.0 );
 106   before->Add(hist); after->Add(hist->Clone()); //13
 107   hist = new TH1D("TPC_vertex_z", "TPC vertex z", 100,-20.0,20.0);
 108   before->Add(hist); after->Add(hist->Clone()); //14
 109   hist = new TH1D("ptyield", "p_{t} spectrum", 10000,0.0,10.0);
 110   before->Add(hist); after->Add(hist->Clone()); //15
 111
 112   //post data here as it doesn't change anyway (the pointer to list anyway)
 113
 114   //restore dir add status
 115   PostData(0, fNtuple);
 116   PostData(1, fOutput);
 117 }
 118
 119 //________________________________________________________________________
 120 void AliAnalysisTaskQAflow::UserExec(Option_t *)
 121 {
 122
 123   //get teh input data
 124   AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
 125   if (!event)
 126   {
 127     AliFatal("no ESD event");
 128     return;
 129   }
 130
 131   //TObjArray* before = ((TDirectory*)fOutput->At(0))->GetList();
 132   //TObjArray* after = ((TDirectory*)fOutput->At(1))->GetList();
 133   TObjArray* before = (TObjArray*)fOutput->At(0);
 134   TObjArray* after = (TObjArray*)fOutput->At(1);
 135   TH1* htrackletmultB = dynamic_cast<TH1*>(before->At(0));
 136   TH1* htrackletmultA = dynamic_cast<TH1*>(after->At(0));
 137   TH1* htrackmultB = dynamic_cast<TH1*>(before->At(1));
 138   TH1* htrackmultA = dynamic_cast<TH1*>(after->At(1));
 139   TH1* hrefmultB = dynamic_cast<TH1*>(before->At(2));
 140   TH1* hrefmultA = dynamic_cast<TH1*>(after->At(2));
 141   TH1* hspdclustersB = dynamic_cast<TH1*>(before->At(3));
 142   TH1* hspdclustersA = dynamic_cast<TH1*>(after->At(3));
 143   TH1* hprimvtxzB = dynamic_cast<TH1*>(before->At(4));
 144   TH1* hprimvtxzA = dynamic_cast<TH1*>(after->At(4));
 145
 146   TH1* hITSclsB = dynamic_cast<TH1*>(before->At(5));
 147   TH1* hITSclsA = dynamic_cast<TH1*>(after->At(5));
 148   TH1* hTPCclsB = dynamic_cast<TH1*>(before->At(6));
 149   TH1* hTPCclsA = dynamic_cast<TH1*>(after->At(6));
 150   TH1* hTPCchi2B = dynamic_cast<TH1*>(before->At(7));
 151   TH1* hTPCchi2A = dynamic_cast<TH1*>(after->At(7));
 152   TH1* hdcaxyB = dynamic_cast<TH1*>(before->At(8));
 153   TH1* hdcaxyA = dynamic_cast<TH1*>(after->At(8));
 154   TH1* hdcazB = dynamic_cast<TH1*>(before->At(9));
 155   TH1* hdcazA = dynamic_cast<TH1*>(after->At(9));
 156   TH1* hphitrackletsB = dynamic_cast<TH1*>(before->At(10));
 157   TH1* hphitrackletsA = dynamic_cast<TH1*>(after->At(10));
 158   TH1* hphitracksB = dynamic_cast<TH1*>(before->At(11));
 159   TH1* hphitracksA = dynamic_cast<TH1*>(after->At(11));
 160   TH1* hetatrackletsB = dynamic_cast<TH1*>(before->At(12));
 161   TH1* hetatrackletsA = dynamic_cast<TH1*>(after->At(12));
 162   TH1* hetatracksB = dynamic_cast<TH1*>(before->At(13));
 163   TH1* hetatracksA = dynamic_cast<TH1*>(after->At(13));
 164   TH1* hprimvtxzTPCB = dynamic_cast<TH1*>(before->At(14));
 165   TH1* hprimvtxzTPCA = dynamic_cast<TH1*>(after->At(14));
 166   TH1* hptyieldB = dynamic_cast<TH1*>(before->At(15));
 167   TH1* hptyieldA = dynamic_cast<TH1*>(after->At(15));
 168
 169   AliMultiplicity* tracklets = const_cast<AliMultiplicity*>(event->GetMultiplicity());
 170   Int_t ntracklets=0;
 171   Int_t nspdclusters=0;
 172   Int_t nspd1clusters=0;
 173   Int_t ntrackletsA=0;
 174   Int_t ntrackletsC=0;
 175   if (tracklets)
 176   {
 177     ntracklets = tracklets->GetNumberOfTracklets();
 178     nspdclusters = tracklets->GetNumberOfITSClusters(0,1);
 179     nspd1clusters = tracklets->GetNumberOfITSClusters(1);
 180     for (Int_t i=0; i<tracklets->GetNumberOfTracklets(); i++)
 181     {
 182       Bool_t pass=fTrackCuts->IsSelected(tracklets,i);
 183       Float_t phi=tracklets->GetPhi(i);
 184       Float_t eta=tracklets->GetEta(i);
 185       hphitrackletsB->Fill(phi); if (pass) hphitrackletsA->Fill(phi);
 186       hetatrackletsB->Fill(eta); if (pass) hetatrackletsA->Fill(eta);
 187       if (eta>0) ntrackletsC++;
 188       else ntrackletsA++;
 189     }
 190
 191   }
 192   //Int_t trackmult = AliESDtrackCuts::GetReferenceMultiplicity(event,kTRUE);
 193   AliFlowEventCuts newcuts;
 194   newcuts.SetRefMultCuts(fTrackCuts);
 195   Int_t trackmult = newcuts.RefMult(event);
 196   Int_t refmult = fEventCuts->RefMult(event);
 197   Bool_t passevent = fEventCuts->IsSelected(event);
 198   htrackletmultB->Fill(ntracklets); if (passevent) htrackletmultA->Fill(ntracklets);
 199   htrackmultB->Fill(trackmult); if (passevent) htrackmultA->Fill( trackmult);
 200   hrefmultB->Fill(refmult); if (passevent) hrefmultA->Fill( refmult);
 201   hspdclustersB->Fill(nspdclusters); if (passevent) hspdclustersA->Fill( nspdclusters);
 202   const AliVertex* vertex = event->GetPrimaryVertex();
 203   Float_t vtxz=0.0;
 204   Float_t vtxx=0.0;
 205   Float_t vtxy=0.0;
 206   if (vertex)
 207   {
 208     vtxz = vertex->GetZ();
 209     vtxx = vertex->GetX();
 210     vtxy = vertex->GetY();
 211     hprimvtxzB->Fill(vtxz); if (passevent) hprimvtxzA->Fill(vtxz);
 212   }
 213   const AliVertex* vertextpc = event->GetPrimaryVertexTPC();
 214   Float_t vtxTPCx=0.0;
 215   Float_t vtxTPCy=0.0;
 216   Float_t vtxTPCz=0.0;
 217   if (vertextpc)
 218   {
 219     vtxTPCx = vertextpc->GetX();
 220     vtxTPCy = vertextpc->GetY();
 221     vtxTPCz = vertextpc->GetZ();
 222     hprimvtxzTPCB->Fill(vtxTPCz); if (passevent) hprimvtxzTPCA->Fill(vtxTPCz);
 223   }
 224   fTrackCuts->SetEvent(event);
 225   Float_t meanpt=0;
 226   Int_t ntracks=fTrackCuts->GetNumberOfInputObjects();
 227   for (Int_t i=0; i<ntracks; i++)
 228   {
 229     TObject* obj = fTrackCuts->GetInputObject(i);
 230     Bool_t pass = fTrackCuts->IsSelected(obj,i);
 231     Float_t dcaxy=0.0;
 232     Float_t dcaz=0.0;
 233     Float_t tpcchi2=0.0;
 234     Float_t tpcchi2percls=0.0;
 235     Int_t ntpccls=0;
 236     Int_t nitscls=0;
 237     Float_t eta=0.0;
 238     Float_t phi=0.0;
 239     Float_t pt=0.0;
 240     AliESDtrack* track = dynamic_cast<AliESDtrack*>(fTrackCuts->GetTrack());
 241     if (track)
 242     {
 243       track->GetImpactParameters(dcaxy,dcaz);
 244       tpcchi2=track->GetTPCchi2();
 245       ntpccls=track->GetTPCNcls();
 246       eta=track->Eta();
 247       phi=track->Phi();
 248       pt=track->Pt();
 249       meanpt+=pt;
 250       tpcchi2percls= (ntpccls==0)?0.0:tpcchi2/ntpccls;
 251       nitscls = track->GetNcls(0);
 252       hITSclsB->Fill(nitscls); if (pass) hITSclsA->Fill( nitscls);
 253       hTPCclsB->Fill(ntpccls); if (pass) hTPCclsA->Fill( ntpccls);
 254       hTPCchi2B->Fill(tpcchi2percls); if (pass) hTPCchi2A->Fill( tpcchi2percls);
 255       hdcaxyB->Fill(dcaxy); if (pass) hdcaxyA->Fill( dcaxy);
 256       hdcazB->Fill(dcaz); if (pass) hdcazA->Fill(dcaz);
 257       hetatracksB->Fill(eta); if (pass) hetatracksA->Fill(eta);
 258       hphitracksB->Fill(phi); if (pass) hphitracksA->Fill(phi);
 259       hptyieldB->Fill(pt); if (pass) hptyieldA->Fill(pt);
 260     }
 261   }
 262   meanpt = meanpt/ntracks;
 263
 264   ///////////////////////////////////////////////////////////////////////
 265   //the ntuple part/////////////
 266   AliFlowEventSimple* flowevent = dynamic_cast<AliFlowEventSimple*>(GetInputData(1));
 267   if (flowevent && fFillNtuple)
 268   {
 269     AliFlowVector qvec = flowevent->GetQ(2);
 270     Double_t qx = qvec.X();
 271     Double_t qy = qvec.Y();
 272     TString triggers = event->GetFiredTriggerClasses();
 273     Int_t trig=0;
 274     if (triggers.Contains("CMBAC-B-NOPF-ALL")) trig+=1;
 275     if (triggers.Contains("CMBS2C-B-NOPF-ALL")) trig+=10;
 276     if (triggers.Contains("CMBS2A-B-NOPF-ALL")) trig+=100;
 277     fNtuple->Fill(meanpt,qx,qy,trackmult,refmult,trig,vtxTPCx,vtxTPCy,vtxTPCz,ntrackletsA,ntrackletsC);
 278   }//if flowevent
 279 }
 280
 281 //________________________________________________________________________
 282 void AliAnalysisTaskQAflow::Terminate(Option_t *)
 283 {
 284
 285 }
 286
 287 //________________________________________________________________________
 288 AliAnalysisTaskQAflow::~AliAnalysisTaskQAflow()
 289 {
 290   //dtor
 291   delete fTrackCuts;
 292   delete fEventCuts;
 293 }