PWG0/FirstPhysics/AliAnalysisTaskFPexample.cxx

   1
   2 #include "AliAnalysisTaskFPexample.h"
   3
   4 #include "TH1D.h"
   5 #include "TH2D.h"
   6 #include "TCanvas.h"
   7 #include "TList.h"
   8 #include "TString.h"
   9
  10 #include "AliAnalysisManager.h"
  11 #include "AliESDtrackCuts.h"
  12 #include "AliESDEvent.h"
  13 #include "AliESDHeader.h"
  14 #include "AliMCEvent.h"
  15 #include "AliMultiplicity.h"
  16 #include "AliTriggerAnalysis.h"
  17
  18 ClassImp(AliAnalysisTaskFPexample)
  19
  20 //________________________________________________________________________
  21 AliAnalysisTaskFPexample::AliAnalysisTaskFPexample(const char *name) :
  22   AliAnalysisTaskFirstPhysics(name),
  23   fhTrackPt(0),
  24   fh2TrackPhiEta(0),
  25   fhMulITSTPC(0),
  26   fhMulITSSA(0),
  27   fhMulSPD(0),
  28   fh2TrackletsPhiEta(0),
  29   fh2TracksPhiTPCchi2(0)
  30 {
  31   DefineOutput(1, TList::Class());
  32 }
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskFPexample::~AliAnalysisTaskFPexample()
  36 {
  37   // Destructor. Clean-up the output list, but not the histograms that are put inside
  38   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
  39   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  40     delete fOutput;
  41   }
  42   for (Int_t i = 0; i < knTrackCuts; i ++) {
  43     delete fTrackCuts[i];
  44     fTrackCuts[i] = 0;
  45   }
  46 }
  47
  48 //________________________________________________________________________
  49 void AliAnalysisTaskFPexample::UserCreateOutputObjects()
  50 {
  51   PrepareOutputList();
  52   // Define cuts
  53   PrepareDefaultTrackCuts();
  54
  55   // Create histograms
  56   const Int_t ptbins = 15;
  57   const Double_t ptlow = 0.1, ptup = 3.1;
  58   const Int_t etabins = 40;
  59   const Double_t etalow = -2.0, etaup = 2.0;
  60   const Int_t mulbins = 200;
  61   const Int_t phibins = 120;
  62
  63   fhTrackPt = UserHisto1d("fhTrackPt", "p_{T} distribution for ESD",
  64                           "P_{T} (GeV/c)", ptbins, ptlow, ptup);
  65   fh2TrackPhiEta = UserHisto2d("fh2TrackPhiEta", "ESD tracks",
  66                                "#Phi", phibins, 0, 2 * TMath::Pi(),
  67                                "#eta", etabins, etalow, etaup);
  68   fhMulITSTPC = UserHisto1d("fhMulITSTPC", "N_{CH} distribution for ITS+TPC tracks",
  69                             "N_{CH}", mulbins, 0, mulbins);
  70   fhMulITSSA = UserHisto1d("fhMulITSSA", "N_{CH} distribution for ITS SA tracks",
  71                            "N_{CH}", mulbins, 0, mulbins);
  72   fhMulSPD = UserHisto1d("fhMulSPD", "N_{CH} distribution for SPD tracklets",
  73                          "N_{CH}", mulbins, 0, mulbins);
  74   fh2TrackletsPhiEta = UserHisto2d("fh2TrackletsPhiEta", "Tracklets",
  75                                    "#Phi", phibins, 0, 2 * TMath::Pi(),
  76                                    "#eta", etabins, etalow, etaup);
  77   fh2TracksPhiTPCchi2 = UserHisto2d("fh2TracksPhiTPCchi2", "ESD tracks #Phi vs. TPC #chi^{2}",
  78                                     "#Phi", phibins, 0, 2 * TMath::Pi(),
  79                                     "#chi^{2}", 500, 0, 500);
  80
  81   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
  82 }
  83
  84 //________________________________________________________________________
  85 void AliAnalysisTaskFPexample::UserExec(Option_t *)
  86 {
  87   if (!GetESDEvent()) {
  88     AliError("Unable to read the ESD");
  89     return;
  90   }
  91   if (!CheckVertex()) {
  92     return;
  93   }
  94
  95   const Int_t nESDTracks = fESD->GetNumberOfTracks();
  96
  97   Int_t iHighestID = 0;
  98   for (Int_t iTrack = 0; iTrack < nESDTracks; iTrack++) {
  99     if (fESD->GetTrack(iTrack)->GetLabel() > iHighestID) {
 100       iHighestID = fESD->GetTrack(iTrack)->GetLabel();
 101     }
 102   }
 103   const Int_t nMaxID = iHighestID + 1;
 104   bool aGlobalBits[nMaxID], aPureITSBits[nMaxID];
 105   for (Int_t iFlag = 0; iFlag < nMaxID; iFlag++) {
 106     aGlobalBits[iFlag] = false;
 107     aPureITSBits[iFlag] = false;
 108   }
 109
 110   // flags for secondary and rejected tracks
 111   const int kRejectBit = BIT(15); // set this bit in ESD tracks if it is rejected by a cut
 112   const int kSecondaryBit = BIT(16); // set this bit in ESD tracks if it is secondary according to a cut
 113
 114   Int_t nTracksITSTPC = 0; // multiplicity counters
 115   Int_t nTracksITSSA = 0;
 116   Int_t nTrackletsSPD = 0;
 117
 118   for(Int_t iTrack = 0; iTrack < nESDTracks; iTrack++){
 119     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
 120     // if track is a secondary from a V0, flag as a secondary
 121     if (esdTrack->IsOn(AliESDtrack::kMultInV0)) {
 122       esdTrack->SetBit(kSecondaryBit);
 123       continue;
 124     }
 125     AliESDtrack* tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD, esdTrack->GetID());
 126     // check tracks with ITS part
 127     if (esdTrack->IsOn(AliESDtrack::kITSin)) {
 128       if (!esdTrack->IsOn(AliESDtrack::kITSpureSA)) { // track has ITS part but is not an ITS_SA
 129         // TPC+ITS
 130         if (esdTrack->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
 131           if (fTrackCuts[kTrackCutQGlo]->AcceptTrack(esdTrack)) { // good ITS+TPC track
 132             if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
 133                 || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
 134               nTracksITSTPC++;
 135               fhTrackPt->Fill(esdTrack->Pt());
 136               fh2TrackPhiEta->Fill(esdTrack->Phi(), esdTrack->Eta());
 137               if (tpcTrack) {
 138                 fh2TracksPhiTPCchi2->Fill(esdTrack->Phi(), tpcTrack->GetTPCchi2());
 139               }
 140             } else {
 141               esdTrack->SetBit(kSecondaryBit); // large DCA -> secondary, don't count either track not associated tracklet
 142             }
 143           } else {
 144             esdTrack->SetBit(kRejectBit); // bad quality, don't count the track, but may count tracklet if associated
 145           }
 146         } else if (fTrackCuts[kTrackCutQITS]->AcceptTrack(esdTrack)) { // good ITS complementary track
 147           if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
 148               || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
 149             nTracksITSTPC++;
 150             fhTrackPt->Fill(esdTrack->Pt());
 151             fh2TrackPhiEta->Fill(esdTrack->Phi(), esdTrack->Eta());
 152             if (tpcTrack) {
 153               fh2TracksPhiTPCchi2->Fill(esdTrack->Phi(), tpcTrack->GetTPCchi2());
 154             }
 155           } else {
 156             esdTrack->SetBit(kSecondaryBit); // large DCA -> secondary, don't count either track not associated tracklet
 157           }
 158         } else {
 159           esdTrack->SetBit(kRejectBit); // bad quality, don't count the track, but may count tracklet if associated
 160         }
 161       } else { // pure ITS SA tracks
 162         if (fTrackCuts[kTrackCutQITS]->AcceptTrack(esdTrack)) { // good ITSSA track
 163           if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
 164               || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
 165             nTracksITSSA++;
 166           } else {
 167             esdTrack->SetBit(kRejectBit);
 168           }
 169         } else {
 170           esdTrack->SetBit(kRejectBit);
 171         }
 172       }
 173     }
 174     if (tpcTrack) {
 175       delete tpcTrack;
 176     }
 177   }
 178
 179   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITS pure SA
 180   const AliMultiplicity* multiplicitySPD = fESD->GetMultiplicity();    // spd multiplicity object
 181   Int_t id1, id2, id3, id4;
 182   AliESDtrack *tr1 = 0, *tr3 = 0;
 183   for (Int_t iTracklet = 0; iTracklet < multiplicitySPD->GetNumberOfTracklets(); iTracklet++) {
 184     if (TMath::Abs(multiplicitySPD->GetEta(iTracklet)) > GetCutEta()) {
 185       continue; // eta selection for tracklets
 186     }
 187     nTrackletsSPD++;
 188     fh2TrackletsPhiEta->Fill(multiplicitySPD->GetPhi(iTracklet), multiplicitySPD->GetEta(iTracklet));
 189     // if counting tracks+tracklets, check if clusters were already used in tracks
 190     // and get the id of the tracks in which they were used
 191
 192     // references for eventual Global/ITS_SA tracks
 193     multiplicitySPD->GetTrackletTrackIDs(iTracklet, 0, id1, id2);
 194     tr1 = id1 >= 0 ? fESD->GetTrack(id1) : 0;
 195
 196     // references for eventual ITS_SA_pure tracks
 197     multiplicitySPD->GetTrackletTrackIDs(iTracklet, 1, id3, id4);
 198     tr3 = id3 >= 0 ? fESD->GetTrack(id3) : 0;
 199
 200     // are both clusters from the same tracks? If not, skip the
 201     // tracklet (shouldn't change things much)
 202     if (id1 != id2 || id3 != id4) {
 203       continue;
 204     }
 205
 206     // has associated global track been associated to a previous tracklet?
 207     bool bUsedInGlobal = (id1 != -1) ? aGlobalBits[id1] : false;
 208     // has associated pure ITS track been associated to a previous tracklet?
 209     bool bUsedInPureITS = (id3 != -1) ? aPureITSBits[id3] : false;
 210
 211     // counting tracklet as global+complementary track
 212     if ((tr1 && !tr1->TestBit(kSecondaryBit)) // reject as secondary
 213         && (tr1 && tr1->TestBit(kRejectBit))) { // already accounted
 214       if (!bUsedInGlobal) {
 215         nTracksITSTPC++;
 216         if (id1 > 0) {
 217           // mark global track linked to this tracklet as associated
 218           aGlobalBits[id1] = true;
 219         }
 220       }
 221     } else if (id1 < 0) {
 222       nTracksITSTPC++;
 223     }
 224
 225     // counting tracklet as ITS SA pure track
 226     if ((tr3 && tr3->TestBit(kSecondaryBit))
 227         && (tr3 && !tr3->TestBit(kRejectBit))) {
 228       if (!bUsedInPureITS) {
 229         nTracksITSSA++;
 230         if (id3 > 0) {
 231           // mark global track linked to this tracklet as associated
 232           aPureITSBits[id3] = true;
 233         }
 234       }
 235     } else if (id3 < 0) {
 236       nTracksITSSA++;
 237     }
 238
 239   }
 240   fhMulITSTPC->Fill(nTracksITSTPC);
 241   fhMulITSSA->Fill(nTracksITSSA);
 242   fhMulSPD->Fill(nTrackletsSPD);
 243
 244   PostData(1, fOutput);
 245 }
 246
 247 //________________________________________________________________________
 248 void AliAnalysisTaskFPexample::Terminate(Option_t *)
 249 {
 250   // Draw result to screen, or perform fitting, normalizations
 251   // don't get too fancy, keep your histos raw and massage them with macros
 252
 253   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 254   if(!fOutput) {
 255     Printf("ERROR: could not retrieve TList fOutput");
 256     return;
 257   }
 258
 259   if (!(GetHisto1FromOutput("fhTrackPt", fhTrackPt) &&
 260         GetHisto2FromOutput("fh2TrackPhiEta", fh2TrackPhiEta) &&
 261         GetHisto1FromOutput("fhMulITSTPC", fhMulITSTPC) &&
 262         GetHisto1FromOutput("fhMulITSSA", fhMulITSSA) &&
 263         GetHisto1FromOutput("fhMulSPD", fhMulSPD) &&
 264         GetHisto2FromOutput("fh2TrackletsPhiEta", fh2TrackletsPhiEta) &&
 265         GetHisto2FromOutput("fh2TracksPhiTPCchi2", fh2TracksPhiTPCchi2))) {
 266     AliError("Couldn't load every histogram from output.");
 267     return;
 268   }
 269
 270   TCanvas *c = new TCanvas("AliAnalysisTaskFPexample", "Data Quality Quick Overview");
 271   c->Divide(2, 2);
 272   c->cd(1)->SetLogy();
 273   fhTrackPt->DrawCopy("E");
 274   c->cd(2);
 275   fh2TrackPhiEta->DrawCopy("");
 276   c->cd(3)->SetLogy();
 277   fhMulITSTPC->DrawCopy("E");
 278   c->cd(4)->SetLogy();
 279   fhMulITSSA->DrawCopy("");
 280 }