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