PWGPP/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 #include "AliESDVZERO.h"
  35 #include "AliMultiplicity.h"
  36 #include "AliPID.h"
  37
  38
  39 #include "AliAnalysisTaskGlobalQA.h"
  40
  41 //
  42 // Run the GlobalQA analysis over already existing ESDs
  43 // Origin: Iouri.Belikov@cern.ch
  44 //
  45
  46 ClassImp(AliAnalysisTaskGlobalQA)
  47
  48 //________________________________________________________________________
  49 AliAnalysisTaskGlobalQA::AliAnalysisTaskGlobalQA() :
  50  AliAnalysisTaskSE("GlobalQA"),
  51  fArrayQA(0)
  52 {
  53   // Default Constructor
  54
  55   DefineOutput(1, TObjArray::Class());
  56 }
  57
  58 //________________________________________________________________________
  59 void AliAnalysisTaskGlobalQA::UserCreateOutputObjects()
  60 {
  61   // Create the histograms
  62   // Called once
  63
  64   fArrayQA = new TObjArray(kLast);
  65
  66
  67   {// Event related QA
  68     const Char_t *name[]={
  69       "hGlobalPrimaryVertex"
  70     };
  71     const Char_t *title[]={
  72       "Z-distribution of the primary vertex"
  73     };
  74     Add2ESDsList(new TH1F(name[0],title[0],100,-20.,20.),kEvt0);
  75   }
  76
  77   {// Cluster related QA
  78     const Char_t *name[]={
  79       "hGlobalFractionAssignedClustersITS",
  80       "hGlobalFractionAssignedClustersTPC",
  81       "hGlobalFractionAssignedClustersTRD",
  82       "hGlobalClustersPerITSModule"
  83      };
  84   const Char_t *title[]={
  85     "Fraction of the assigned clusters in ITS",
  86     "Fraction of the assigned clusters in TPC",
  87     "Fraction of the assigned clusters in TRD",
  88     "Number of clusters per an ITS module"
  89   };
  90   Add2ESDsList(new TH1F(name[0],title[0],100,0.,2.),kClr0);
  91   Add2ESDsList(new TH1F(name[1],title[1],100,0.,2.),kClr1);
  92   Add2ESDsList(new TH1F(name[2],title[2],100,0.,2.),kClr2);
  93   Add2ESDsList(new TH1F(name[3],title[3],2201,-0.5,2200.5),kClr3);
  94   }
  95
  96   {// Track related QA
  97     const Char_t *name[]={
  98       "hGlobalTrackAzimuthe",                               // kTrk0
  99       "hGlobalTrackEta",                                    // kTrk1
 100       "hGlobalTPCTrackpT",                                  // kTrk2
 101       "hGlobalTPCITSMatchedpT",                             // kTrk3
 102       "hGlobalTPCTOFMatchedpT",                             // kTrk4
 103       "hGlobalTPCITSMatchingProbability",                   // kTrk5
 104       "hGlobalTPCTOFMatchingProbability",                   // kTrk6
 105       "hGlobalTPCsideAposDCA",                              // kTrk7
 106       "hGlobalTPCsideAnegDCA",                              // kTrk8
 107       "hGlobalTPCsideCposDCA",                              // kTrk9
 108       "hGlobalTPCsideCnegDCA"                               // kTrk10
 109   };
 110   const Char_t *title[]={
 111     "Track azimuthal distribution (rad)",                   // kTrk0
 112     "Track pseudo-rapidity distribution",                   // kTrk1
 113     "TPC: track momentum distribution (GeV)",               // kTrk2
 114     "TPC-ITS matched: track momentum distribution (GeV)",   // kTrk3
 115     "TPC-TOF matched: track momentum distribution (GeV)",   // kTrk4
 116     "TPC-ITS track-matching probability",                   // kTrk5
 117     "TPC-TOF track-matching probability",                   // kTrk6
 118     "TPC side A: DCA for the positive tracks (mm)",         // kTrk7
 119     "TPC side A: DCA for the negative tracks (mm)",         // kTrk8
 120     "TPC side C: DCA for the positive tracks (mm)",         // kTrk9
 121     "TPC side C: DCA for the negative tracks (mm)"          // kTrk10
 122   };
 123   Add2ESDsList(new TH1F(name[0],title[0],100, 0.,TMath::TwoPi()),kTrk0);
 124   Add2ESDsList(new TH1F(name[1],title[1],100,-2.00,2.00),kTrk1);
 125   Add2ESDsList(new TH1F(name[2],title[2],50,  0.20,5.00),kTrk2);
 126   Add2ESDsList(new TH1F(name[3],title[3],50,  0.20,5.00),kTrk3);
 127   Add2ESDsList(new TH1F(name[4],title[4],50,  0.20,5.00),kTrk4);
 128   Add2ESDsList(new TH1F(name[5],title[5],50,  0.20,5.00),kTrk5);
 129   Add2ESDsList(new TH1F(name[6],title[6],50,  0.20,5.00),kTrk6);
 130   Add2ESDsList(new TH1F(name[7],title[7],50, -25.0,25.0),kTrk7);
 131   Add2ESDsList(new TH1F(name[8],title[8],50, -25.0,25.0),kTrk8);
 132   Add2ESDsList(new TH1F(name[9],title[9],50, -25.0,25.0),kTrk9);
 133   Add2ESDsList(new TH1F(name[10],title[10],50, -25.0,25.0),kTrk10);
 134   }
 135
 136   {// V0 related QA
 137     const Char_t *name[]={
 138       "hGlobalPromptK0sMass",
 139       "hGlobalOfflineK0sMass",
 140       "hGlobalPromptLambda0Lambda0BarMass",
 141       "hGlobalOfflineLambda0Lambda0BarMass"
 142     };
 143   const Char_t *title[]={
 144     "On-the-fly K0s mass (GeV)",
 145     "Offline K0s mass (GeV)",
 146     "On-the-fly Lambda0 + Lambda0Bar mass (GeV)",
 147     "Offline Lambda0 + Lambda0Bar mass (GeV)"
 148   };
 149   Add2ESDsList(new TH1F(name[0],title[0],50,  0.4477,0.5477),kK0on);
 150   Add2ESDsList(new TH1F(name[1],title[1],50,  0.4477,0.5477),kK0off);
 151   Add2ESDsList(new TH1F(name[2],title[2],50,  1.0657,1.1657),kL0on);
 152   Add2ESDsList(new TH1F(name[3],title[3],50,  1.0657,1.1657),kL0off);
 153   }
 154
 155   {// PID related QA
 156   const Char_t *name[]={
 157     "hGlobalITSdEdx",
 158     "hGlobalTPCdEdx",
 159     "hGlobalTOFTrackingvsMeasured",
 160     "hGlobalTPCdEdxvsMomentum"
 161   };
 162   const Char_t *title[]={
 163     "ITS: dEdx (A.U.) for particles with momentum 0.4 - 0.5 (GeV)",
 164     "TPC: dEdx (A.U.) for particles with momentum 0.4 - 0.5 (GeV)",
 165     "TOF: tracking - measured (ps)",
 166     "TPC: dEdx (A.U.) vs momentum (GeV)"
 167   };
 168   Add2ESDsList(new TH1F(name[0],title[0],50,0.00,200.),kPid0);
 169   Add2ESDsList(new TH1F(name[1],title[1],50,0.00,100.),kPid1);
 170   Add2ESDsList(new TH1F(name[2],title[2],50,-3500.,3500.),kPid2);
 171   Add2ESDsList(new TH2F(name[3],title[3],1500,0.05,15.,700,0.,700.),kPid3);
 172   }
 173
 174   {// Multiplicity related QA
 175     const Char_t *name[]={
 176       "hGlobalV0AvsITS",
 177       "hGlobalV0CvsITS"
 178     };
 179     const Char_t *title[]={
 180       "Multiplicity: V0A vs ITS",
 181       "Multiplicity: V0C vs ITS"
 182     };
 183     TH2F *h0=new TH2F(name[0],title[0],41,-0.5,40.5, 33,-0.5,32.5);
 184     Add2ESDsList(h0,kMlt0);
 185     TH2F *h1=new TH2F(name[1],title[1],41,-0.5,40.5, 33,-0.5,32.5);
 186     Add2ESDsList(h1,kMlt1);
 187   }
 188
 189 }
 190
 191 //________________________________________________________________________
 192 void AliAnalysisTaskGlobalQA::UserExec(Option_t *)
 193 {
 194   // Main loop
 195   // Called for each event
 196
 197   const AliESDEvent *esd=(const AliESDEvent *)InputEvent();
 198
 199   if (!esd) {
 200     Printf("ERROR: ESD is not available");
 201     return;
 202   }
 203
 204   // Event related QA
 205   const AliESDVertex *vtx=esd->GetPrimaryVertex();
 206   if (!vtx->GetStatus()) return;
 207
 208   Double_t xv=vtx->GetXv();
 209   Double_t yv=vtx->GetYv();
 210   Double_t zv=vtx->GetZv();
 211   GetESDsData(kEvt0)->Fill(zv);
 212
 213
 214   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 215       AliESDtrack* track = esd->GetTrack(iTracks);
 216       if (!track) {
 217          Printf("ERROR: Could not receive track %d", iTracks);
 218          continue;
 219       }
 220
 221
 222     // Cluster related QA
 223     if (track->IsOn(AliESDtrack::kITSrefit)) {
 224       Int_t n=track->GetITSclusters(0);
 225       GetESDsData(kClr0)->Fill(Float_t(n)/6.); //6 is the number of ITS layers
 226     }
 227
 228     for (Int_t i=0; i<6; i++) {
 229       Int_t idet, sts;
 230       Float_t xloc,zloc;
 231       if (!track->GetITSModuleIndexInfo(i,idet,sts,xloc,zloc)) continue;
 232       if (i>=2) idet+=240;
 233       if (i>=4) idet+=260;
 234       if ((sts==1)||(sts==2)||(sts==4)) GetESDsData(kClr3)->Fill(idet);
 235     }
 236
 237     if (track->IsOn(AliESDtrack::kTPCrefit)) {
 238       Int_t n =track->GetTPCNcls();
 239       Int_t nf=track->GetTPCNclsF();      // number of crossed TPC pad rows
 240       if (nf>0) {
 241         Double_t val = n*1.0/nf;
 242         GetESDsData(kClr1)->Fill(val);
 243       }
 244     }
 245
 246     if (track->IsOn(AliESDtrack::kTRDrefit)) {
 247       Int_t n=track->GetTRDclusters(0);
 248       GetESDsData(kClr2)->Fill(Float_t(n)/(6*24));//(6*24) is the number of TRD time bins
 249     }
 250
 251     Double_t p=track->GetP();
 252
 253     // Track related QA
 254     if (track->IsOn(AliESDtrack::kTPCrefit)) {
 255       Float_t dz[2];
 256       track->GetDZ(xv,yv,zv,esd->GetMagneticField(),dz);
 257       if ((TMath::Abs(dz[0])<3.) && (TMath::Abs(dz[1])<3.)) { // beam pipe
 258         Double_t phi=track->Phi();
 259         GetESDsData(kTrk0)->Fill(phi);
 260         Double_t y=track->Eta();
 261         GetESDsData(kTrk1)->Fill(y);
 262
 263         if (TMath::Abs(y)<0.9) {
 264            GetESDsData(kTrk2)->Fill(p);
 265            if (track->IsOn(AliESDtrack::kITSrefit)) GetESDsData(kTrk3)->Fill(p);
 266           //if (track->IsOn(AliESDtrack::kTOFout)) GetESDsData(kTrk4)->Fill(p);
 267            if (track->GetTOFsignal()>0) GetESDsData(kTrk4)->Fill(p);
 268         }
 269       }
 270     }
 271     const AliExternalTrackParam *tpcTrack=track->GetTPCInnerParam();
 272     const AliExternalTrackParam *innTrack=track->GetInnerParam();
 273     if (tpcTrack)
 274     if (innTrack) {
 275        Float_t dz[2];
 276        tpcTrack->GetDZ(xv,yv,zv,esd->GetMagneticField(),dz);
 277        dz[0]*=10.; // in mm
 278        if (innTrack->GetZ()  > 0)
 279        if (innTrack->GetTgl()> 0) { // TPC side A
 280           if (tpcTrack->GetSign() > 0) GetESDsData(kTrk7)->Fill(dz[0]);
 281           else                         GetESDsData(kTrk8)->Fill(dz[0]);
 282        }
 283        if (innTrack->GetZ()  < 0)
 284        if (innTrack->GetTgl()< 0) { // TPC side C
 285           if (tpcTrack->GetSign() > 0) GetESDsData(kTrk9)->Fill(dz[0]);
 286           else                         GetESDsData(kTrk10)->Fill(dz[0]);
 287        }
 288     }
 289
 290     // PID related QA
 291     if ((p>0.4) && (p<0.5)) {
 292       if (track->IsOn(AliESDtrack::kITSpid)) {
 293         Double_t dedx=track->GetITSsignal();
 294         GetESDsData(kPid0)->Fill(dedx);
 295       }
 296       if (track->IsOn(AliESDtrack::kTPCpid)) {
 297         Double_t dedx=track->GetTPCsignal();
 298         GetESDsData(kPid1)->Fill(dedx);
 299       }
 300     }
 301     if (p>1.0) {
 302       if (track->IsOn(AliESDtrack::kITSrefit))
 303       if (track->IsOn(AliESDtrack::kTPCrefit))
 304       if (track->IsOn(AliESDtrack::kTOFout)) {
 305          Float_t dz[2];
 306          track->GetDZ(xv,yv,zv,esd->GetMagneticField(),dz);
 307          if (dz[1]<3.) {
 308            Double_t times[AliPID::kSPECIESC];
 309            track->GetIntegratedTimes(times,AliPID::kSPECIESC);
 310            Double_t tof=track->GetTOFsignal()/*-847055 -1771207*/;
 311            GetESDsData(kPid2)->Fill(times[2]-tof);
 312          }
 313       }
 314     }
 315     const AliExternalTrackParam *par=track->GetInnerParam();
 316     if (par) {
 317       Double_t pp=par->GetP();
 318       Double_t dedx=track->GetTPCsignal();
 319       TH2F *h = dynamic_cast<TH2F*>(GetESDsData(kPid3));
 320       if (h) h->Fill(pp,dedx);
 321     }
 322   }
 323
 324   // Multiplicity related QA
 325   AliESDVZERO     *mltV0 =esd->GetVZEROData();
 326   const AliMultiplicity *mltITS=esd->GetMultiplicity();
 327   if (mltV0)
 328     if (mltITS) {
 329        Short_t nv0a=mltV0->GetNbPMV0A();
 330        Short_t nv0c=mltV0->GetNbPMV0C();
 331        Int_t   nits=mltITS->GetNumberOfTracklets();
 332        TH2F *h0=dynamic_cast<TH2F*>(GetESDsData(kMlt0));
 333        if (h0) h0->Fill(nits,nv0a);
 334        TH2F *h1=dynamic_cast<TH2F*>(GetESDsData(kMlt1));
 335        if (h1) h1->Fill(nits,nv0c);
 336     }
 337
 338   // V0 related QA
 339   Int_t nV0=esd->GetNumberOfV0s();
 340   for (Int_t i=0; i<nV0; i++) {
 341     Double_t mass;
 342     AliESDv0 v0(*esd->GetV0(i));
 343
 344     Int_t nidx=TMath::Abs(v0.GetNindex());
 345     AliESDtrack *ntrack1=esd->GetTrack(nidx);
 346     if (!ntrack1->IsOn(AliESDtrack::kTPCrefit)) continue;
 347
 348     Int_t pidx=TMath::Abs(v0.GetPindex());
 349     AliESDtrack *ptrack1=esd->GetTrack(pidx);
 350     if (!ptrack1->IsOn(AliESDtrack::kTPCrefit)) continue;
 351
 352     v0.ChangeMassHypothesis(kK0Short);
 353     mass=v0.GetEffMass();
 354     if (v0.GetOnFlyStatus())
 355        GetESDsData(kK0on)->Fill(mass);
 356     else
 357        GetESDsData(kK0off)->Fill(mass);
 358
 359     v0.ChangeMassHypothesis(kLambda0);
 360     mass=v0.GetEffMass();
 361     if (v0.GetOnFlyStatus())
 362        GetESDsData(kL0on)->Fill(mass);
 363     else
 364        GetESDsData(kL0off)->Fill(mass);
 365
 366     v0.ChangeMassHypothesis(kLambda0Bar);
 367     mass=v0.GetEffMass();
 368     if (v0.GetOnFlyStatus())
 369        GetESDsData(kL0on)->Fill(mass);
 370     else
 371        GetESDsData(kL0off)->Fill(mass);
 372   }
 373
 374   // Post output data.
 375   PostData(1, fArrayQA);
 376 }
 377
 378 //________________________________________________________________________
 379 void AliAnalysisTaskGlobalQA::Terminate(Option_t *)
 380 {
 381   // Draw the results on the screen
 382   // Called once at the end of the query
 383
 384   fArrayQA=(TObjArray*)GetOutputData(1);
 385
 386   TH1 *tpc=GetESDsData(kTrk2); tpc->Sumw2();
 387   TH1 *its=GetESDsData(kTrk3); its->Sumw2();
 388   TH1 *tof=GetESDsData(kTrk4); tof->Sumw2();
 389   GetESDsData(kTrk5)->Divide(its,tpc,1,1.,"b");
 390   GetESDsData(kTrk6)->Divide(tof,tpc,1,1.,"b");
 391
 392   TH1 *hTPCdEdxMIP = GetESDsData(kPid1);
 393   if (!hTPCdEdxMIP) {
 394     Printf("ERROR: hTPCdEdxMIP not available");
 395     return;
 396   }
 397
 398   TH2 *hTPCdEdxVsP = dynamic_cast<TH2*>(GetESDsData(kPid3));
 399   if (!hTPCdEdxVsP) {
 400     Printf("ERROR: hTPCdEdxVsP not available");
 401     return;
 402   }
 403
 404   TCanvas *c2=new TCanvas("c2","",320,32,530,590);
 405
 406   TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw(); p6->cd();
 407   p6->SetFillColor(42); p6->SetFrameFillColor(10);
 408   hTPCdEdxMIP->SetFillColor(2); hTPCdEdxMIP->SetFillStyle(3005);
 409   if (hTPCdEdxMIP->GetEntries()<333)
 410       hTPCdEdxMIP->DrawCopy("E");
 411   else
 412       hTPCdEdxMIP->Fit("gaus");
 413   c2->cd();
 414
 415   TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw(); p7->cd(); p7->SetLogx();
 416   p7->SetFillColor(42); p7->SetFrameFillColor(10);
 417   hTPCdEdxVsP->DrawCopy();
 418
 419 }