PWG1/global/AliAnalysisTaskGlobalQA.cxx

   1
   2 /**************************************************************************
   3  *  Authors : Iouri Belikov, Antonin Maire
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //-----------------------------------------------------------------
  16 //                 AliAnalysisTaskGlobalQA class
  17 // This task is for running the GlobalQA over already existing ESDs
  18 //          Origin:  I.Belikov, Iouri.Belikov@cern.ch, June 2009
  19 //-----------------------------------------------------------------
  20
  21 #include <TPDGCode.h>
  22
  23 #include "TChain.h"
  24 #include "TTree.h"
  25 #include "TH1F.h"
  26 #include "TH2F.h"
  27 #include "TCanvas.h"
  28
  29 #include "AliAnalysisManager.h"
  30
  31 #include "AliESDEvent.h"
  32 #include "AliESDv0.h"
  33 #include "AliESDInputHandler.h"
  34
  35 #include "AliAnalysisTaskGlobalQA.h"
  36
  37 //
  38 // Run the GlobalQA analysis over already existing ESDs
  39 // Origin: Iouri.Belikov@cern.ch
  40 //
  41
  42 ClassImp(AliAnalysisTaskGlobalQA)
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskGlobalQA::AliAnalysisTaskGlobalQA() :
  46  AliAnalysisTaskSE("GlobalQA"),
  47  fArrayQA(0)
  48 {
  49   // Default Constructor
  50
  51   DefineOutput(1, TObjArray::Class());
  52 }
  53
  54 //________________________________________________________________________
  55 void AliAnalysisTaskGlobalQA::UserCreateOutputObjects()
  56 {
  57   // Create the histograms
  58   // Called once
  59
  60   fArrayQA = new TObjArray(kLast);
  61
  62
  63   {// Cluster related QA
  64     const Char_t *name[]={
  65       "hGlobalFractionAssignedClustersITS",
  66       "hGlobalFractionAssignedClustersTPC",
  67       "hGlobalFractionAssignedClustersTRD"
  68     };
  69   const Char_t *title[]={
  70     "Fraction of the assigned clusters in ITS",
  71     "Fraction of the assigned clusters in TPC",
  72     "Fraction of the assigned clusters in TRD"
  73   };
  74   Add2ESDsList(new TH1F(name[0],title[0],100,0.,2.),kClr0);
  75   Add2ESDsList(new TH1F(name[1],title[1],100,0.,2.),kClr1);
  76   Add2ESDsList(new TH1F(name[2],title[2],100,0.,2.),kClr2);
  77   }
  78
  79   {// Track related QA
  80     const Char_t *name[]={
  81       "hGlobalTrackAzimuthe",                               // kTrk0
  82       "hGlobalTrackEta",                                    // kTrk1
  83       "hGlobalTPCTrackpT",                                  // kTrk2
  84       "hGlobalTPCITSMatchedpT",                             // kTrk3
  85       "hGlobalTPCTOFMatchedpT",                             // kTrk4
  86       "hGlobalTPCITSMatchingProbability",                   // kTrk5
  87       "hGlobalTPCTOFMatchingProbability",                   // kTrk6
  88       "hGlobalTPCsideAposDCA",                              // kTrk7
  89       "hGlobalTPCsideAnegDCA",                              // kTrk8
  90       "hGlobalTPCsideCposDCA",                              // kTrk9
  91       "hGlobalTPCsideCnegDCA"                               // kTrk10
  92   };
  93   const Char_t *title[]={
  94     "Track azimuthal distribution (rad)",                   // kTrk0
  95     "Track pseudo-rapidity distribution",                   // kTrk1
  96     "TPC: track momentum distribution (GeV)",               // kTrk2
  97     "TPC-ITS matched: track momentum distribution (GeV)",   // kTrk3
  98     "TPC-TOF matched: track momentum distribution (GeV)",   // kTrk4
  99     "TPC-ITS track-matching probability",                   // kTrk5
 100     "TPC-TOF track-matching probability"                    // kTrk6
 101     "TPC side A: DCA for the positive tracks (mm)",         // kTrk7
 102     "TPC side A: DCA for the negative tracks (mm)",         // kTrk8
 103     "TPC side C: DCA for the positive tracks (mm)",         // kTrk9
 104     "TPC side C: DCA for the negative tracks (mm)"          // kTrk10
 105   };
 106   Add2ESDsList(new TH1F(name[0],title[0],100, 0.,TMath::TwoPi()),kTrk0);
 107   Add2ESDsList(new TH1F(name[1],title[1],100,-2.00,2.00),kTrk1);
 108   Add2ESDsList(new TH1F(name[2],title[2],50,  0.20,5.00),kTrk2);
 109   Add2ESDsList(new TH1F(name[3],title[3],50,  0.20,5.00),kTrk3);
 110   Add2ESDsList(new TH1F(name[4],title[4],50,  0.20,5.00),kTrk4);
 111   Add2ESDsList(new TH1F(name[5],title[5],50,  0.20,5.00),kTrk5);
 112   Add2ESDsList(new TH1F(name[6],title[6],50,  0.20,5.00),kTrk6);
 113   Add2ESDsList(new TH1F(name[7],title[7],50, -25.0,25.0),kTrk7);
 114   Add2ESDsList(new TH1F(name[8],title[8],50, -25.0,25.0),kTrk8);
 115   Add2ESDsList(new TH1F(name[9],title[9],50, -25.0,25.0),kTrk9);
 116   Add2ESDsList(new TH1F(name[10],title[10],50, -25.0,25.0),kTrk10);
 117   }
 118
 119   {// V0 related QA
 120     const Char_t *name[]={
 121       "hGlobalPromptK0sMass",
 122       "hGlobalOfflineK0sMass",
 123       "hGlobalPromptLambda0Lambda0BarMass",
 124       "hGlobalOfflineLambda0Lambda0BarMass"
 125     };
 126   const Char_t *title[]={
 127     "On-the-fly K0s mass (GeV)",
 128     "Offline K0s mass (GeV)",
 129     "On-the-fly Lambda0 + Lambda0Bar mass (GeV)",
 130     "Offline Lambda0 + Lambda0Bar mass (GeV)"
 131   };
 132   Add2ESDsList(new TH1F(name[0],title[0],50,  0.4477,0.5477),kK0on);
 133   Add2ESDsList(new TH1F(name[1],title[1],50,  0.4477,0.5477),kK0off);
 134   Add2ESDsList(new TH1F(name[2],title[2],50,  1.0657,1.1657),kL0on);
 135   Add2ESDsList(new TH1F(name[3],title[3],50,  1.0657,1.1657),kL0off);
 136   }
 137
 138   {// PID related QA
 139   const Char_t *name[]={
 140     "hGlobalITSdEdx",
 141     "hGlobalTPCdEdx",
 142     "hGlobalTOFTrackingvsMeasured",
 143     "hGlobalTPCdEdxvsMomentum"
 144   };
 145   const Char_t *title[]={
 146     "ITS: dEdx (A.U.) for particles with momentum 0.4 - 0.5 (GeV)",
 147     "TPC: dEdx (A.U.) for particles with momentum 0.4 - 0.5 (GeV)",
 148     "TOF: tracking - measured (ps)",
 149     "TPC: dEdx (A.U.) vs momentum (GeV)"
 150   };
 151   Add2ESDsList(new TH1F(name[0],title[0],50,0.00,200.),kPid0);
 152   Add2ESDsList(new TH1F(name[1],title[1],50,0.00,100.),kPid1);
 153   Add2ESDsList(new TH1F(name[2],title[2],50,-3500.,3500.),kPid2);
 154   Add2ESDsList(new TH2F(name[3],title[3],1500,0.05,15.,700,0.,700.),kPid3);
 155   }
 156
 157 }
 158
 159 //________________________________________________________________________
 160 void AliAnalysisTaskGlobalQA::UserExec(Option_t *)
 161 {
 162   // Main loop
 163   // Called for each event
 164
 165   const AliESDEvent *esd=(const AliESDEvent *)InputEvent();
 166
 167   if (!esd) {
 168     Printf("ERROR: ESD is not available");
 169     return;
 170   }
 171
 172   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 173       AliESDtrack* track = esd->GetTrack(iTracks);
 174       if (!track) {
 175          Printf("ERROR: Could not receive track %d", iTracks);
 176          continue;
 177       }
 178
 179
 180     // Cluster related QA
 181     if (track->IsOn(AliESDtrack::kITSrefit)) {
 182       Int_t n=track->GetITSclusters(0);
 183       GetESDsData(kClr0)->Fill(Float_t(n)/6.); //6 is the number of ITS layers
 184     }
 185
 186     if (track->IsOn(AliESDtrack::kTPCrefit)) {
 187       Int_t n =track->GetTPCNcls();
 188       Int_t nf=track->GetTPCNclsF();      // number of crossed TPC pad rows
 189       if (nf>0) {
 190         Double_t val = n*1.0/nf;
 191         GetESDsData(kClr1)->Fill(val);
 192       }
 193     }
 194
 195     if (track->IsOn(AliESDtrack::kTRDrefit)) {
 196       Int_t n=track->GetTRDclusters(0);
 197       GetESDsData(kClr2)->Fill(Float_t(n)/(6*24));//(6*24) is the number of TRD time bins
 198     }
 199
 200     Double_t p=track->GetP();
 201
 202     // Track related QA
 203     if (track->IsOn(AliESDtrack::kTPCrefit)) {
 204       Float_t dz[2];
 205       track->GetDZ(0.,0.,0.,esd->GetMagneticField(),dz);
 206       if ((TMath::Abs(dz[0])<3.) && (TMath::Abs(dz[1])<3.)) { // beam pipe
 207         Double_t phi=track->Phi();
 208         GetESDsData(kTrk0)->Fill(phi);
 209         Double_t y=track->Eta();
 210         GetESDsData(kTrk1)->Fill(y);
 211
 212         if (TMath::Abs(y)<0.9) {
 213            GetESDsData(kTrk2)->Fill(p);
 214            if (track->IsOn(AliESDtrack::kITSrefit)) GetESDsData(kTrk3)->Fill(p);
 215           //if (track->IsOn(AliESDtrack::kTOFout)) GetESDsData(kTrk4)->Fill(p);
 216            if (track->GetTOFsignal()>0) GetESDsData(kTrk4)->Fill(p);
 217         }
 218       }
 219     }
 220     const AliExternalTrackParam *tpcTrack=track->GetTPCInnerParam();
 221     const AliExternalTrackParam *innTrack=track->GetInnerParam();
 222     if (tpcTrack)
 223     if (innTrack) {
 224        const AliESDVertex *vtx=esd->GetPrimaryVertex();
 225        Double_t xv=vtx->GetXv();
 226        Double_t yv=vtx->GetYv();
 227        Double_t zv=vtx->GetZv();
 228        Float_t dz[2];
 229        tpcTrack->GetDZ(xv,yv,zv,esd->GetMagneticField(),dz);
 230        dz[0]*=10.; // in mm
 231        if (innTrack->GetZ()  > 0)
 232        if (innTrack->GetTgl()> 0) { // TPC side A
 233           if (tpcTrack->GetSign() > 0) GetESDsData(kTrk7)->Fill(dz[0]);
 234           else                         GetESDsData(kTrk8)->Fill(dz[0]);
 235        }
 236        if (innTrack->GetZ()  < 0)
 237        if (innTrack->GetTgl()< 0) { // TPC side C
 238           if (tpcTrack->GetSign() > 0) GetESDsData(kTrk9)->Fill(dz[0]);
 239           else                         GetESDsData(kTrk10)->Fill(dz[0]);
 240        }
 241     }
 242
 243     // PID related QA
 244     if ((p>0.4) && (p<0.5)) {
 245       if (track->IsOn(AliESDtrack::kITSpid)) {
 246         Double_t dedx=track->GetITSsignal();
 247         GetESDsData(kPid0)->Fill(dedx);
 248       }
 249       if (track->IsOn(AliESDtrack::kTPCpid)) {
 250         Double_t dedx=track->GetTPCsignal();
 251         GetESDsData(kPid1)->Fill(dedx);
 252       }
 253     }
 254     if (p>1.0) {
 255       if (track->IsOn(AliESDtrack::kTOFpid)) {
 256         Double_t times[10];
 257         track->GetIntegratedTimes(times);
 258         Double_t tof=track->GetTOFsignal();
 259         GetESDsData(kPid2)->Fill(times[2]-tof);
 260       }
 261     }
 262     const AliExternalTrackParam *par=track->GetInnerParam();
 263     if (par) {
 264       Double_t pp=par->GetP();
 265       Double_t dedx=track->GetTPCsignal();
 266       TH2F *h = dynamic_cast<TH2F*>(GetESDsData(kPid3));
 267       h->Fill(pp,dedx);
 268     }
 269   }
 270
 271   // V0 related QA
 272   Int_t nV0=esd->GetNumberOfV0s();
 273   for (Int_t i=0; i<nV0; i++) {
 274     Double_t mass;
 275     AliESDv0 v0(*esd->GetV0(i));
 276
 277     v0.ChangeMassHypothesis(kK0Short);
 278     mass=v0.GetEffMass();
 279     if (v0.GetOnFlyStatus())
 280        GetESDsData(kK0on)->Fill(mass);
 281     else
 282        GetESDsData(kK0off)->Fill(mass);
 283
 284     v0.ChangeMassHypothesis(kLambda0);
 285     mass=v0.GetEffMass();
 286     if (v0.GetOnFlyStatus())
 287        GetESDsData(kL0on)->Fill(mass);
 288     else
 289        GetESDsData(kL0off)->Fill(mass);
 290
 291     v0.ChangeMassHypothesis(kLambda0Bar);
 292     mass=v0.GetEffMass();
 293     if (v0.GetOnFlyStatus())
 294        GetESDsData(kL0on)->Fill(mass);
 295     else
 296        GetESDsData(kL0off)->Fill(mass);
 297   }
 298
 299   // Post output data.
 300   PostData(1, fArrayQA);
 301 }
 302
 303 //________________________________________________________________________
 304 void AliAnalysisTaskGlobalQA::Terminate(Option_t *)
 305 {
 306   // Draw the results on the screen
 307   // Called once at the end of the query
 308
 309   fArrayQA=(TObjArray*)GetOutputData(1);
 310
 311   TH1 *tpc=GetESDsData(kTrk2); tpc->Sumw2();
 312   TH1 *its=GetESDsData(kTrk3); its->Sumw2();
 313   TH1 *tof=GetESDsData(kTrk4); tof->Sumw2();
 314   GetESDsData(kTrk5)->Divide(its,tpc,1,1.,"b");
 315   GetESDsData(kTrk6)->Divide(tof,tpc,1,1.,"b");
 316
 317   TH1 *hTPCdEdxMIP = GetESDsData(kPid1);
 318   if (!hTPCdEdxMIP) {
 319     Printf("ERROR: hTPCdEdxMIP not available");
 320     return;
 321   }
 322
 323   TH2 *hTPCdEdxVsP = dynamic_cast<TH2*>(GetESDsData(kPid3));
 324   if (!hTPCdEdxVsP) {
 325     Printf("ERROR: hTPCdEdxVsP not available");
 326     return;
 327   }
 328
 329   TCanvas *c2=new TCanvas("c2","",320,32,530,590);
 330
 331   TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw(); p6->cd();
 332   p6->SetFillColor(42); p6->SetFrameFillColor(10);
 333   hTPCdEdxMIP->SetFillColor(2); hTPCdEdxMIP->SetFillStyle(3005);
 334   if (hTPCdEdxMIP->GetEntries()<333)
 335       hTPCdEdxMIP->DrawCopy("E");
 336   else
 337       hTPCdEdxMIP->Fit("gaus");
 338   c2->cd();
 339
 340   TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw(); p7->cd(); p7->SetLogx();
 341   p7->SetFillColor(42); p7->SetFrameFillColor(10);
 342   hTPCdEdxVsP->DrawCopy();
 343
 344 }