TPC/Calib/AliAnalysisTaskPt.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TCanvas.h"
   5 #include "TObjArray.h"
   6
   7 #include "AliAnalysisTask.h"
   8 #include "AliAnalysisManager.h"
   9
  10 #include "AliESDEvent.h"
  11 #include "AliESDtrackCuts.h"
  12 #include "AliVEventHandler.h"
  13 #include "AliTPCseed.h"
  14 #include "AliTPCclusterMI.h"
  15 #include "AliVfriendEvent.h"
  16 #include "AliVfriendTrack.h"
  17 #include "AliESDInputHandler.h"
  18
  19 #include "AliAnalysisTaskPt.h"
  20
  21 // example of an analysis task creating a p_t spectrum
  22 // Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing
  23
  24 ClassImp(AliAnalysisTaskPt)
  25
  26 //________________________________________________________________________
  27 AliAnalysisTaskPt::AliAnalysisTaskPt(const char *name)
  28 : AliAnalysisTask(name, ""), fESD(0), fESDfriend(0), fHistPt(0), fCuts(0), fEv(0), fHistQ(0), fListOut(0), fUseFriends(kFALSE), fHistNTPCCl(0), fHistNESDtracks(0),   fHistNESDfriendtracks(0)
  29
  30 {
  31   // Constructor
  32
  33   // Define input and output slots here
  34   // Input slot #0 works with a TChain
  35   DefineInput(0, TChain::Class());
  36   // Output slot #0 writes into a TH1 container
  37   DefineOutput(0, TList::Class());
  38 }
  39
  40 //________________________________________________________________________
  41 void AliAnalysisTaskPt::ConnectInputData(Option_t *)
  42 {
  43   // Connect ESD or AOD here
  44   // Called once
  45
  46   printf("----> AliAnalysisTaskPt::ConnectInputData\n");
  47   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  48   if (!tree) {
  49     Printf("ERROR: Could not read chain from input slot 0");
  50   }
  51   else {
  52     // Disable all branches and enable only the needed ones
  53     // The next two lines are different when data produced as AliESDEvent is read
  54     /*
  55     tree->SetBranchStatus("*", kFALSE);
  56     tree->SetBranchStatus("fTracks.*", kTRUE);
  57     */
  58
  59     AliVEventHandler *esdH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  60     TString classInputHandler = esdH->ClassName();
  61
  62     Printf("----> AliAnalysisTaskPt: ClassName of handler = %s", classInputHandler.Data());
  63
  64     if (!esdH) {
  65       Printf("ERROR: Could not get ESDInputHandler");
  66     }
  67     else {
  68       Printf("----> AliAnalysisTaskPt::ConnectInputData Getting the event from handler %p", esdH);
  69       fESD = esdH->GetEvent();
  70       if (fUseFriends){
  71         Printf("...We have to use the friends...");
  72         if (classInputHandler.Contains("HLT")) { // we are running in HLT
  73           fESDfriend = esdH->GetVFriendEvent();
  74         }
  75         else { /// we are running offline
  76           if (esdH && esdH->GetTree()) {
  77             Printf("...We got the tree...");
  78             if (esdH->GetTree()->GetBranch("ESDfriend.")){
  79               Printf("Yu-huuuu!!! friend branch found");
  80               fESDfriend = ((AliESDInputHandler*)esdH)->GetESDfriend();
  81             }
  82             else {
  83               Printf("No friend branch found");
  84             }
  85           }
  86         }
  87
  88         Printf("and the result is: fESDfriend = %p", fESDfriend);
  89       }
  90       else {
  91         Printf("The friends are not requested");
  92       }
  93     }
  94     if (!fESD) {
  95       Printf("ERROR, no ESD event");
  96     }
  97     if (fUseFriends && !fESDfriend) {
  98       Printf("ERROR, no ESD friend");
  99     }
 100   }
 101
 102   Printf("fESD = %p, fESDfriend = %p", fESD, fESDfriend);
 103   printf("<---- AliAnalysisTaskPt::ConnectInputData\n");
 104 }
 105
 106 //________________________________________________________________________
 107 void AliAnalysisTaskPt::CreateOutputObjects()
 108 {
 109   // Create histograms
 110   // Called once
 111
 112   fListOut = new TList();
 113   fListOut->SetOwner();
 114   fListOut->SetName("listHistos");
 115
 116   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
 117   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 118   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 119   fHistPt->SetMarkerStyle(kFullCircle);
 120
 121   fHistQ = new TH1F("fHistQ", "TPC clusters: Q distribution", 1000, 0, 10000);
 122   fHistQ->GetXaxis()->SetTitle("Q");
 123   fHistQ->GetYaxis()->SetTitle("dN/dQ");
 124   fHistQ->SetMarkerStyle(kFullCircle);
 125
 126   fHistNTPCCl = new TH1F("fHistNTPCCl", "Number of TPC clusters", 160, -0.5, 159.5);
 127   fHistNTPCCl->GetXaxis()->SetTitle("n. TPC Cl.");
 128   fHistNTPCCl->GetYaxis()->SetTitle("dN/d(n. TPC Cl)");
 129   fHistNTPCCl->SetMarkerStyle(kFullCircle);
 130
 131   fHistNESDtracks = new TH1F("fHistNESDtracks", "Number of ESD friend tracks", 1000, -0.5, 999.5);
 132   fHistNESDtracks->GetXaxis()->SetTitle("n. ESD friend tracks");
 133   fHistNESDtracks->GetYaxis()->SetTitle("dN/d(n. ESD friend tracks)");
 134   fHistNESDtracks->SetMarkerStyle(kFullCircle);
 135
 136   fHistNESDfriendtracks = new TH1F("fHistNESDfriendtracks", "Number of ESD tracks", 1000, -0.5, 999.5);
 137   fHistNESDfriendtracks->GetXaxis()->SetTitle("n. ESD tracks");
 138   fHistNESDfriendtracks->GetYaxis()->SetTitle("dN/d(n. ESD tracks)");
 139   fHistNESDfriendtracks->SetMarkerStyle(kFullCircle);
 140
 141   fListOut->Add(fHistPt);
 142   fListOut->Add(fHistQ);
 143   fListOut->Add(fHistNTPCCl);
 144   fListOut->Add(fHistNESDtracks);
 145   fListOut->Add(fHistNESDfriendtracks);
 146
 147   PostData(0, fListOut);
 148
 149   fCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
 150 }
 151
 152 //________________________________________________________________________
 153 void AliAnalysisTaskPt::Exec(Option_t *)
 154 {
 155   // Main loop
 156   // Called for each event
 157
 158   if (!fESD) {
 159     Printf("ERROR: fESD not available");
 160     return;
 161   }
 162
 163   /*
 164   if (fUseFriends){
 165     Printf("In Exec: ...We have to use the friends...");
 166     fESDfriend = fESD->FindFriend();
 167     Printf("...and we got friends = %p", fESDfriend);
 168     if (!fESDfriend) {
 169       Printf("ERROR: fESDfriend not available");
 170         return;
 171     }
 172   }
 173   */
 174
 175   if (fUseFriends){
 176     Printf("In Exec: ...We have to use the friends...");
 177     Printf("...and we got friends = %p", fESDfriend);
 178     if (!fESDfriend) {
 179       Printf("ERROR: fESDfriend not available");
 180         return;
 181     }
 182   }
 183
 184   Int_t nESDtracks = fESD->GetNumberOfTracks();
 185   Int_t nESDfriendtracks = 0;
 186   if (fUseFriends) nESDfriendtracks = fESDfriend->GetNumberOfTracks();
 187   Printf("There are %d tracks in this event", nESDtracks);
 188   Printf("... and there are %d friends in this event", nESDfriendtracks);
 189
 190   fHistNESDtracks->Fill(nESDtracks);
 191   fHistNESDfriendtracks->Fill(nESDfriendtracks);
 192
 193   // Track loop to fill a pT spectrum
 194   for (Int_t iTracks = 0; iTracks < nESDtracks; iTracks++) {
 195     Printf("Checking track %d", iTracks);
 196     const AliVTrack* track = dynamic_cast<AliVTrack*>(fESD->GetTrack(iTracks));
 197     if (!track) {
 198       Printf("ERROR: Could not receive track %d", iTracks);
 199       continue;
 200     }
 201     Printf("track %d has pt = %f", iTracks, track->Pt());
 202     fHistPt->Fill(track->Pt());
 203     fHistNTPCCl->Fill(track->GetTPCNcls());
 204   } //track loop
 205
 206
 207   if (fUseFriends){
 208     Printf("In the loop over the friends");
 209     // Friend Track loop
 210     for (Int_t iFriend = 0; iFriend < nESDfriendtracks; iFriend++) {
 211       Printf("Getting friend %d", iFriend);
 212       const AliVfriendTrack* friendTrack = fESDfriend->GetTrack(iFriend);
 213       if (!friendTrack) {
 214         Printf("ERROR: Could not receive track %d", iFriend);
 215         continue;
 216       }
 217       else {
 218         Printf("friend track = %p", friendTrack);
 219       }
 220
 221       AliTPCseed seed;
 222       Int_t err = friendTrack->GetTPCseed( seed );
 223       Printf("err = %d", err);
 224       if( err==0 ){
 225         Printf("Found TPC seed" );
 226         for (Int_t irow = 0; irow < 160; irow++){
 227           AliTPCclusterMI* cluMI = seed.GetClusterPointer(irow);
 228           if (cluMI){
 229             Printf("Found cluster at row %d", irow);
 230             Float_t q = cluMI->GetQ();
 231             Printf("Q = %f", q);
 232             fHistQ->Fill(q);
 233           }
 234           else {
 235             Printf("Row %d does not contain clusters", irow);
 236           }
 237         }
 238       }
 239       else {
 240         //Printf("Schade... seed is %p", seed);
 241       }
 242     }
 243   }
 244
 245   // Post output data.
 246   PostData(0, fListOut);
 247   fEv++;
 248 }
 249
 250 //________________________________________________________________________
 251 void AliAnalysisTaskPt::Terminate(Option_t *)
 252 {
 253   // Draw result to the screen
 254   // Called once at the end of the query
 255
 256   Printf("Terminate called: fESD = %p", fESD);
 257
 258   fListOut = dynamic_cast<TList*> (GetOutputData(0));
 259
 260   if (fListOut) {
 261     fHistPt = dynamic_cast<TH1F*>(fListOut->FindObject("fHistPt"));
 262     if (!fHistPt) {
 263       Printf("ERROR: fHistPt not available");
 264       return;
 265     }
 266
 267     TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
 268     c1->cd(1)->SetLogy();
 269     fHistPt->DrawCopy("E");
 270   }
 271   else {
 272     Printf("In Terminate: no TList found");
 273   }
 274
 275 }