PWGLF/SPECTRA/ChargedHadrons/multPbPb/AliAnalysisTaskTriggerStudy.cxx

   1 // AliAnalysisTaskTriggerStudy
   2
   3 // Author: Michele Floris, CERN
   4 // TODO:
   5 // - Add chi2/cluster plot for primary, secondaries and fakes
   6
   7
   8 #include "AliAnalysisTaskTriggerStudy.h"
   9 #include "AliESDInputHandler.h"
  10 #include "AliHistoListWrapper.h"
  11 #include "AliAnalysisManager.h"
  12 #include "AliMCEvent.h"
  13 #include "AliStack.h"
  14 #include "TH1I.h"
  15 #include "TH3D.h"
  16 #include "AliMCParticle.h"
  17 #include "AliGenEventHeader.h"
  18 #include "AliCentrality.h"
  19
  20 #include <iostream>
  21 #include "AliTriggerAnalysis.h"
  22 #include "AliMultiplicity.h"
  23 #include "TFile.h"
  24 #include "AliLog.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliESDVZERO.h"
  27 #include "TH2F.h"
  28 #include "AliESDUtils.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliGenDPMjetEventHeader.h"
  31
  32 using namespace std;
  33
  34 ClassImp(AliAnalysisTaskTriggerStudy)
  35
  36 //const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"C0SM1","C0SM2","C0VBA","C0VBC","C0OM2"};
  37 const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"V0AND online","V0AND offline","Physics Selection", "Rec Candle", "Gen Canlde"};//,"C0OM2"};
  38
  39 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy()
  40 : AliAnalysisTaskSE("TaskTriggerStudy"),
  41   fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  42 {
  43   // constructor
  44
  45   DefineOutput(1, AliHistoListWrapper::Class());
  46
  47 }
  48 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const char * name)
  49   : AliAnalysisTaskSE(name),
  50     fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  51 {
  52   //
  53   // Standard constructur which should be used
  54   //
  55
  56   DefineOutput(1, AliHistoListWrapper::Class());
  57
  58 }
  59
  60 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const AliAnalysisTaskTriggerStudy& obj) :
  61   AliAnalysisTaskSE(obj) ,fESD (0), fIsMC(0), fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  62 {
  63   //copy ctor
  64   fESD = obj.fESD ;
  65   fHistoList = obj.fHistoList;
  66   fTriggerAnalysis = obj.fTriggerAnalysis;
  67   fHistoSuffix = obj.fHistoSuffix;
  68   fNTrackletsCut = obj.fNTrackletsCut;
  69   fNTrackletsCutKine = obj.fNTrackletsCutKine;
  70   fRejectBGWithV0 = obj.fRejectBGWithV0;
  71 }
  72
  73 AliAnalysisTaskTriggerStudy::~AliAnalysisTaskTriggerStudy(){
  74   // destructor
  75
  76   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  77     if(fHistoList) {
  78       delete fHistoList;
  79       fHistoList = 0;
  80     }
  81     if(fTriggerAnalysis) {
  82       delete fTriggerAnalysis;
  83       fHistoList = 0;
  84     }
  85   }
  86   // Histo list should not be destroyed: fListWrapper is owner!
  87
  88 }
  89 void AliAnalysisTaskTriggerStudy::UserCreateOutputObjects()
  90 {
  91   // Called once
  92   fHistoList = new AliHistoListWrapper("histoList","histogram list for trigger studies");
  93   fTriggerAnalysis = new AliTriggerAnalysis();
  94   if (fIsMC) fTriggerAnalysis->SetAnalyzeMC(1);
  95 }
  96
  97
  98 void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
  99 {
 100   // User code
 101
 102   // FIXME: make sure you have the right cuts here
 103
 104   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 105   PostData(1,fHistoList);
 106
 107   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 108   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
 109     AliFatal("Not processing ESDs");
 110   }
 111
 112   // get the multiplicity object
 113   const AliMultiplicity* mult = fESD->GetMultiplicity();
 114   Int_t ntracklets = mult->GetNumberOfTracklets();
 115   // Get Number of tracks
 116   Int_t ntracks    = AliESDtrackCuts::GetReferenceMultiplicity(fESD,kTRUE); // tpc only
 117
 118   // Get V0 Multiplicity
 119   AliESDVZERO* esdV0 = fESD->GetVZEROData();
 120   Float_t multV0A=esdV0->GetMTotV0A();
 121   Float_t multV0C=esdV0->GetMTotV0C();
 122   Float_t dummy = 0;
 123   Float_t multV0 =  AliESDUtils::GetCorrV0(fESD, dummy);
 124
 125   // Get number of clusters in layer 1
 126   Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);
 127   Float_t innerLayerSPD = mult->GetNumberOfITSClusters(0);
 128   Float_t totalClusSPD = outerLayerSPD+innerLayerSPD;
 129
 130   //  if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CMBS2A-B-NOPF-ALL")|| fESD->IsTriggerClassFired("CMBS2C-B-NOPF-ALL") || fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL")) )) return;
 131   //  if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")))) return;
 132   if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))) return;
 133   GetHistoSPD1  ("All", "All events before any selection")->Fill(outerLayerSPD);
 134   GetHistoV0M   ("All", "All events before any selection")->Fill(multV0);
 135
 136   // V0 triggers
 137   Bool_t c0v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
 138   Bool_t c0v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
 139   // Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 140   // Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 141   Bool_t v0AHW = fTriggerAnalysis->EvaluateTrigger(fESD, AliTriggerAnalysis::kCTPV0A);
 142   Bool_t v0CHW = fTriggerAnalysis->EvaluateTrigger(fESD, AliTriggerAnalysis::kCTPV0C);
 143   Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0ABG);
 144   Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0CBG);
 145   Bool_t v0BG  = v0ABG || v0CBG;
 146
 147   // At least one track in eta < 0.8 with pt > 0.5
 148   // MC Checks
 149   Bool_t atLeast1Track = kFALSE;
 150   Bool_t isSD = kFALSE;
 151   Bool_t isND = kFALSE;
 152
 153   if(fIsMC) {
 154     if (!fMCEvent) {
 155       AliError("No MC info found");
 156     } else {
 157       Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
 158       Int_t nPhysicalPrimaries = 0;
 159       for (Int_t ipart=0; ipart<nMCTracks; ipart++) {
 160         AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 161         // We don't care about neutrals and non-physical primaries
 162         if(mcPart->Charge() == 0) continue;
 163
 164         //check if current particle is a physical primary
 165         if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
 166         if(mcPart->Pt()<0.5) continue;
 167         if(TMath::Abs(mcPart->Eta())>0.8) continue;
 168         atLeast1Track = kTRUE;
 169         break;
 170       }
 171
 172       AliGenPythiaEventHeader * headPy  = 0;
 173       AliGenDPMjetEventHeader * headPho = 0;
 174       AliGenEventHeader * htmp = fMCEvent->GenEventHeader();
 175       if(!htmp) {
 176         AliFatal("Cannot Get MC Header!!");
 177         return;
 178       }
 179       if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
 180         headPy =  (AliGenPythiaEventHeader*) htmp;
 181       } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
 182         headPho = (AliGenDPMjetEventHeader*) htmp;
 183       } else {
 184         AliFatal("Unknown header");
 185       }
 186       if(headPy)   {
 187         //    cout << "Process: " << headPy->ProcessType() << endl;
 188         if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
 189           isSD = kTRUE; // is single difractive
 190         }
 191         if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
 192           isND = kTRUE; // is non-diffractive
 193         }
 194
 195       } else if (headPho) {
 196         if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
 197           isSD = kTRUE;
 198         }
 199         if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
 200           isND = kTRUE;
 201         }
 202       }
 203
 204     }
 205   }
 206
 207   static AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 208   // cuts->SetPtRange(0.5,10000);
 209   // cuts->SetEtaRange(-0.8, 0.8);
 210   Int_t ntracksLoop = fESD->GetNumberOfTracks();
 211
 212   Bool_t atLeast1TrackESD = kFALSE;
 213   for (Int_t iTrack = 0; iTrack<ntracksLoop; iTrack++) {
 214     AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
 215     // for each track
 216
 217     // track quality cuts
 218     if(!cuts->AcceptTrack(track)) continue;
 219     if(track->Pt()<0.5) continue;
 220     if(TMath::Abs(track->Eta())>0.8) continue;
 221
 222     atLeast1TrackESD = kTRUE;
 223     break;
 224   }
 225
 226
 227
 228   // Physics selection
 229   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 230
 231   Bool_t vdArray[kNVDEntries];
 232   vdArray[kVDV0ANDOnline]      = v0AHW && v0CHW && !v0BG ;
 233   vdArray[kVDV0ANDOffline]     = c0v0A && c0v0C && isSelected;
 234   vdArray[kVDRecCandle]        = atLeast1TrackESD;
 235   vdArray[kVDPhysSel]          = isSelected;
 236   vdArray[kVDGenCandle]        = atLeast1Track;
 237
 238   FillTriggerOverlaps("All", "All Events",vdArray);
 239   if(!isSD)   FillTriggerOverlaps("NSD", "NSD Events",vdArray);
 240
 241
 242   if(!isSelected) return;
 243
 244
 245   GetHistoSPD1  ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(outerLayerSPD);
 246   GetHistoV0M   ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(multV0);
 247
 248
 249   // Francesco's cuts
 250   // const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
 251   // if(vtxESDTPC->GetNContributors()<1) return;
 252   // if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;
 253   // const AliESDVertex * vtxESDSPD= fESD->GetPrimaryVertexSPD();
 254   // Float_t tpcContr=vtxESDTPC->GetNContributors();
 255
 256
 257
 258   // GetT0 Stuff
 259   const Double32_t *meanT0   = fESD->GetT0TOF();
 260   const Double32_t  meanT0A  = 0.001* meanT0[1];
 261   const Double32_t  meanT0C  = 0.001* meanT0[2];
 262   const Double32_t  meanT0AC = 0.001* meanT0[0];
 263   Double32_t  T0Vertex = fESD->GetT0zVertex();
 264   //  Double32_t  *ampT0 =Esdevent ->GetT0amplitude();
 265
 266 //   cut1yesF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && francescoscut
 267 // cut1notF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && ! francescoscut
 268 // cut2 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) <=-0.7) && meanC > -2) )
 269 // cut3 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) < 0.7 && (meanC-meanA) > -0.7 ) )
 270 //   cut4 = ( (meanC < 95. && meanA < 95.) && (meanA < -2.)
 271
 272   Bool_t cut1T0 =  ( (meanT0C < 95. && meanT0A < 95.) && (meanT0C < -2.) );
 273   Bool_t cut2T0 = ( (meanT0C < 95. && meanT0A < 95.) && ( (meanT0C-meanT0A) <=-0.7) && meanT0C > -2) ;
 274
 275   if(ntracklets > fNTrackletsCut) return;
 276
 277   // Reset histo suffix and fill reference histograms without any suffix
 278   fHistoSuffix = "";
 279
 280   // Fast or in the outer layer
 281   Int_t nFastOrOnline  = fTriggerAnalysis->SPDFiredChips(fESD, 1, 0, 2); // offline
 282   Int_t nFastOrOffline = fTriggerAnalysis->SPDFiredChips(fESD, 0, 0, 2); // online
 283
 284   Bool_t c0sm1 = nFastOrOffline >= 1;
 285   Bool_t c0sm2 = nFastOrOffline >= 2;
 286   Bool_t c0sm3 = nFastOrOffline >= 3;
 287   Bool_t c0sm4 = nFastOrOffline >= 4;
 288   Bool_t c0sm5 = nFastOrOffline >= 5;
 289
 290
 291   // TOF triggers
 292   // FIXME: move to triggeranalysis?
 293   AliESDHeader*h = fESD->GetHeader(); // taken the header from AliESDEvent
 294   Bool_t c0OM2 = h->IsTriggerInputFired("0OM2"); // thr >= 2 (input 19)
 295   Bool_t c0OM3 = h->IsTriggerInputFired("0OM3"); // thr >= 3 (input 20)
 296
 297   // ZDC triggers
 298   Bool_t zdcA    = kFALSE;
 299   Bool_t zdcC    = kFALSE;
 300   Bool_t zdcBar  = kFALSE;
 301   Bool_t zdcTime = fTriggerAnalysis->ZDCTimeTrigger(fESD);
 302
 303   if (!fIsMC) {
 304     // If it's data, we use the TDCs
 305     zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ;
 306     zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;
 307     zdcBar = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCentralBarrel) ;
 308   } else {
 309     // If it's MC, we use the energy
 310     Double_t minEnergy = 0;
 311     AliESDZDC *esdZDC = fESD->GetESDZDC();
 312     Double_t  zNCEnergy = esdZDC->GetZDCN1Energy();
 313     Double_t  zPCEnergy = esdZDC->GetZDCP1Energy();
 314     Double_t  zNAEnergy = esdZDC->GetZDCN2Energy();
 315     Double_t  zPAEnergy = esdZDC->GetZDCP2Energy();
 316     // zdcA = (zNAEnergy>minEnergy || zPAEnergy>minEnergy);
 317     // zdcC = (zNCEnergy>minEnergy || zPCEnergy>minEnergy);
 318     zdcA = (zNAEnergy>minEnergy);
 319     zdcC = (zNCEnergy>minEnergy);
 320   }
 321
 322
 323   // Some macros for the online triggers
 324   Bool_t cMBS2A = c0sm2 && c0v0A;
 325   Bool_t cMBS2C = c0sm2 && c0v0C;
 326   Bool_t cMBAC  = c0v0A && c0v0C;
 327
 328
 329
 330
 331
 332   // Plots requested by Jurgen on 18/11/2010 + later updates (including plots for the note)
 333   // FIXME: will skip everything else
 334
 335
 336   if(zdcTime) {
 337     GetHistoSPD1  ("ZDCTIME", "ZDC Time Cut")->Fill(outerLayerSPD);
 338     GetHistoTracks("ZDCTIME", "ZDC Time Cut")->Fill(ntracks);
 339     GetHistoV0M   ("ZDCTIME", "ZDC Time Cut")->Fill(multV0);
 340   }
 341
 342   if(zdcTime && ntracks > 1) {
 343     GetHistoSPD1  ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(outerLayerSPD);
 344     GetHistoTracks("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(ntracks);
 345     GetHistoV0M   ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(multV0);
 346   }
 347
 348   // GetHistoSPD1  ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(outerLayerSPD);
 349   // GetHistoTracks("PhysSel", "All events after physics selection and Francesco's cut")->Fill(ntracks);
 350   // GetHistoV0M   ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(multV0);
 351   if(c0v0A && c0v0C) {
 352     GetHistoSPD1  ("V0AND", "V0A & V0C")->Fill(outerLayerSPD);
 353     GetHistoTracks("V0AND", "V0A & V0C")->Fill(ntracks);
 354     GetHistoV0M   ("V0AND", "V0A & V0C")->Fill(multV0);
 355   }
 356   if((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) {
 357     GetHistoSPD1  ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(outerLayerSPD);
 358     GetHistoTracks("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(ntracks);
 359     GetHistoV0M   ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(multV0);
 360   }
 361   if(zdcA && zdcC) {
 362     GetHistoSPD1  ("ZDCAND", "ZDCA & ZDCC")->Fill(outerLayerSPD);
 363     GetHistoTracks("ZDCAND", "ZDCA & ZDCC")->Fill(ntracks);
 364     GetHistoV0M   ("ZDCAND", "ZDCA & ZDCC")->Fill(multV0);
 365   }
 366   if((c0v0A && c0v0C) && !(zdcA && zdcC)) {
 367     GetHistoSPD1  ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(outerLayerSPD);
 368     GetHistoTracks("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(ntracks);
 369     GetHistoV0M   ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(multV0);
 370   }
 371   if((c0v0A && c0v0C) && !zdcA && !zdcC) {
 372     GetHistoSPD1  ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(outerLayerSPD);
 373     GetHistoTracks("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(ntracks);
 374     GetHistoV0M   ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(multV0);
 375 }
 376   if((c0v0A && c0v0C) && ((!zdcA && zdcC) || (zdcA && !zdcC))) {
 377     GetHistoSPD1  ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(outerLayerSPD);
 378     GetHistoTracks("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(ntracks);
 379     GetHistoV0M   ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(multV0);
 380   }
 381
 382   if(((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) && (zdcA && zdcC)) {
 383     GetHistoSPD1  ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(outerLayerSPD);
 384     GetHistoTracks("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(ntracks);
 385     GetHistoV0M   ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(multV0);
 386   }
 387
 388   if((c0v0A && c0v0C) && (zdcA && zdcC)) {
 389     GetHistoSPD1  ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(outerLayerSPD);
 390     GetHistoTracks("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(ntracks);
 391     GetHistoV0M   ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(multV0);
 392   }
 393
 394   if((c0v0A && zdcA && !zdcC && !c0v0C) || (c0v0C && zdcC && !zdcA && !c0v0A)) {
 395     GetHistoSPD1  ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(outerLayerSPD);
 396     GetHistoTracks("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(ntracks);
 397     GetHistoV0M   ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(multV0);
 398   }
 399
 400   // GetHistoCorrelationSPDTPCVz("All", "After physics selection and Francesco's cut")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 401   // if(cut1T0)   GetHistoCorrelationSPDTPCVz("Cut1T0", "T0 Cut 1")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 402   // if(cut2T0)   GetHistoCorrelationSPDTPCVz("Cut2T0", "T0 Cut 2")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 403
 404   // GetHistoCorrelationContrTPCSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,tpcContr);
 405   // if(cut1T0)   GetHistoCorrelationContrTPCSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,tpcContr);
 406   // if(cut2T0)   GetHistoCorrelationContrTPCSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,tpcContr);
 407
 408   // GetHistoCorrelationTrackletsSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,ntracklets);
 409   // if(cut1T0)   GetHistoCorrelationTrackletsSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,ntracklets);
 410   // if(cut2T0)   GetHistoCorrelationTrackletsSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,ntracklets);
 411
 412
 413   return; // FIXME
 414
 415
 416   // Reject background
 417   if (v0BG && fRejectBGWithV0) {
 418     cout << "Rejection BG" << endl;
 419
 420     return;
 421   }
 422   // Fill global histos
 423   GetHistoTracklets("all","All events")->Fill(ntracklets);
 424   FillTriggerOverlaps("All", "All Events",vdArray);
 425
 426   // Fill some combination of trigger classes
 427   Bool_t cmbs1aOnline = fESD->IsTriggerClassFired("CMBS1A-B-NOPF-ALL");
 428   Bool_t cmbs1cOnline = fESD->IsTriggerClassFired("CMBS1C-B-NOPF-ALL");
 429   Bool_t cmbacOnline  = fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL");
 430
 431   Bool_t twoOutOfThree = kFALSE;
 432   if (cmbs1aOnline || cmbs1cOnline ||cmbacOnline) twoOutOfThree = kTRUE;
 433   if (twoOutOfThree) GetHistoTracklets("All_TwoOutOfThree" ,"Events 2-out-of-3 online" )->Fill(ntracklets);
 434
 435
 436   // loop over trigger classes in the event
 437   TObjArray * tokens = 0;
 438   if(fIsMC) {
 439     // in case of montecarlo I override the trigger class
 440     tokens = new TObjArray;
 441     tokens->SetOwner();
 442     //    tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL"));
 443     tokens->Add(new TObjString("MC"));
 444   }
 445   else {
 446     TString trgClasses = fESD->GetFiredTriggerClasses();
 447     tokens = trgClasses.Tokenize("  ");
 448   }
 449   TIterator  * iter = (TIterator*) tokens->MakeIterator();
 450
 451   TString classes = fESD->GetFiredTriggerClasses();
 452
 453   // if (classes.Contains("SMH")) {
 454   //    tokens->Print();
 455   //    cout << classes.Data() << endl;
 456   // }
 457   //  iter->Reset();
 458   //Int_t itoken = 0;
 459   TObjString * tok=0;
 460   while((tok = (TObjString*) iter->Next())){
 461     // clean up trigger name
 462     TString trg = tok->GetString();
 463     trg.Strip(TString::kTrailing, ' ');
 464     trg.Strip(TString::kLeading, ' ');
 465
 466     fHistoSuffix = "_";
 467     fHistoSuffix += trg;
 468
 469     //    cout << itoken++ << " " << trg.Data() << endl;
 470     // continue;
 471     //    if (trg.Contains("SMH")) cout << itoken++ << " " << trg.Data() << endl;
 472
 473     // Fill tracklets
 474     GetHistoTracklets("All" ,"Events no offline trigger" )->Fill(ntracklets);
 475
 476
 477     // Fill histograms mismatchs
 478     // TODO: check mismatch trigger class
 479     if(nFastOrOffline != nFastOrOnline) {
 480       GetHistoTracklets("mismatchingFastOr", "Events where fast or offline differs from fast-or online")->Fill(ntracklets);
 481     }
 482
 483     if (c0v0A != v0AHW){
 484       GetHistoTracklets("mismatchingV0A", "Events where V0A offline differs from V0A online")->Fill(ntracklets);
 485     }
 486
 487     if (c0v0C != v0CHW){
 488       GetHistoTracklets("mismatchingV0C", "Events where V0C offline differs from V0C online")->Fill(ntracklets);
 489     }
 490
 491     // Fill a tracklet histo for each trigger type
 492     if(c0sm1)  GetHistoTracklets("c0sm1" ,"Events were trigger c0sm1 fired" )->Fill(ntracklets);
 493     if(c0sm2)  GetHistoTracklets("c0sm2" ,"Events were trigger c0sm2 fired" )->Fill(ntracklets);
 494     if(c0sm3)  GetHistoTracklets("c0sm3" ,"Events were trigger c0sm3 fired" )->Fill(ntracklets);
 495     if(c0sm4)  GetHistoTracklets("c0sm4" ,"Events were trigger c0sm4 fired" )->Fill(ntracklets);
 496     if(c0sm5)  GetHistoTracklets("c0sm5" ,"Events were trigger c0sm5 fired" )->Fill(ntracklets);
 497     if(c0OM2)  GetHistoTracklets("c0OM2" ,"Events were trigger c0OM2 fired" )->Fill(ntracklets);
 498     if(c0OM3)  GetHistoTracklets("c0OM3" ,"Events were trigger c0OM3 fired" )->Fill(ntracklets);
 499     if(c0v0A)  GetHistoTracklets("c0v0A" ,"Events were trigger c0v0A fired" )->Fill(ntracklets);
 500     if(c0v0C)  GetHistoTracklets("c0v0C" ,"Events were trigger c0v0C fired" )->Fill(ntracklets);
 501     if(cMBS2A) GetHistoTracklets("cMBS2A","Events were trigger cMBS2A fired")->Fill(ntracklets);
 502     if(cMBS2C) GetHistoTracklets("cMBS2C","Events were trigger cMBS2C fired")->Fill(ntracklets);
 503     if(cMBAC ) GetHistoTracklets("cMBAC", "Events were trigger cMBAC  fired")->Fill(ntracklets);
 504     if(zdcA )  GetHistoTracklets("cZDCA", "Events were trigger cZDCA  fired")->Fill(ntracklets);
 505     if(zdcC )  GetHistoTracklets("cZDCC", "Events were trigger cZDCC  fired")->Fill(ntracklets);
 506     if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC", "Events were zdc trigger don't  fired")->Fill(ntracklets);
 507     if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC", "Events were only 1 ZDC  fired")->Fill(ntracklets);
 508     if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC", "Events were both  ZDC  fired")->Fill(ntracklets);
 509
 510     if(twoOutOfThree) {
 511       if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC_TwoOutOfThree", "Events were zdc trigger don't  fired")->Fill(ntracklets);
 512       if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC_TwoOutOfThree", "Events were only 1 ZDC  fired")->Fill(ntracklets);
 513       if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC_TwoOutOfThree", "Events were both  ZDC  fired")->Fill(ntracklets);
 514     }
 515     if(zdcBar) GetHistoTracklets("cZDCBar","Events were trigger cZDCB  fired")->Fill(ntracklets);
 516     //  if() GetHistoTracklets("","Events were trigger  fired");
 517
 518     // Fill trigger overlaps
 519     FillTriggerOverlaps("All", "All Events in trigger class",vdArray);
 520
 521
 522     // Fill kinematic variables for peripheral events
 523     static AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); // FIXME: change this ?
 524     if (ntracklets < fNTrackletsCutKine) {
 525       // 1. Loop on tracklets
 526       for(Int_t itracklet = 0; itracklet < ntracklets; itracklet++){
 527         GetHistoEta("All", Form("Tracklet #eta distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(mult->GetEta(itracklet));
 528       }
 529       // 2. loop on tracks
 530       TObjArray * goodTracks = cuts->GetAcceptedTracks(fESD);
 531       TIterator * iterTracks = goodTracks->MakeIterator();
 532       AliESDtrack * esdTrack = 0;
 533       while ((esdTrack = (AliESDtrack*) iterTracks->Next())) {
 534         GetHistoPt("All", Form("Tracklet p_{T} distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(esdTrack->Pt());
 535       }
 536       delete goodTracks;
 537     }
 538
 539
 540   }
 541   delete tokens;
 542
 543     // if (fIsMC) {
 544
 545
 546   //   if (!fMCEvent) {
 547   //     AliError("No MC info found");
 548   //   } else {
 549
 550   //     //loop on the MC event
 551   //     //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
 552   //     Int_t offset    = fMCEvent->GetPrimaryOffset();
 553   //     Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
 554   //     for (Int_t ipart=offset; ipart<nMCTracks; ipart++) {
 555
 556   //    AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 557
 558   //    // We don't care about neutrals and non-physical primaries
 559   //    if(mcPart->Charge() == 0) continue;
 560
 561   // PHYSICAL PRIMARY
 562   //    // Get MC vertex
 563     //  TArrayF   vertex;
 564   //    fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
 565   //    Float_t zv = vertex[2];
 566   //    //      Float_t zv = vtxESD->GetZ();
 567   //    // Fill generated histo
 568   //    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
 569
 570   //     }
 571   //   }
 572   // }
 573
 574
 575
 576
 577 }
 578
 579 void   AliAnalysisTaskTriggerStudy::Terminate(Option_t *){
 580   // terminate
 581   // Save output in a more friendly format
 582   fHistoList = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
 583   if (!fHistoList){
 584     Printf("ERROR: fHistoList not available");
 585     return;
 586   }
 587   TFile * f = new TFile("trigger_study.root", "recreate");
 588   fHistoList->GetList()->Write();
 589   f->Close();
 590
 591 }
 592
 593 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracklets(const char * name, const char * title){
 594   // Book histo of events vs ntracklets, if needed
 595
 596   TString hname = "hTracklets_";
 597   hname+=name;
 598   hname+=fHistoSuffix;
 599   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 600
 601   if(!h) {
 602     AliInfo(Form("Booking histo %s",hname.Data()));
 603     Bool_t oldStatus = TH1::AddDirectoryStatus();
 604     TH1::AddDirectory(kFALSE);
 605     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 606     h = new TH1F (hname.Data(), title, 4000, -0.5, 4000-0.5);
 607     h->Sumw2();
 608     h->SetXTitle("ntracklets");
 609     fHistoList->GetList()->Add(h);
 610     TH1::AddDirectory(oldStatus);
 611   }
 612   return h;
 613 }
 614
 615
 616 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoSPD1(const char * name, const char * title){
 617   // Book histo of events vs ntracklets, if needed
 618
 619   TString hname = "hSPD1_";
 620   hname+=name;
 621   hname+=fHistoSuffix;
 622   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 623
 624   // const Int_t nbin = 118;
 625   // Float_t bins[nbin] = {0};
 626   // bins[0] = 0;
 627   // for(Int_t ibin = 1; ibin <= nbin; ibin++){
 628   //   if (ibin < 100)
 629   //     bins[ibin] = bins [ibin-1]+1;
 630   //   else if (ibin < 1000)
 631   //     bins[ibin] = bins [ibin-1]+100;
 632   //   else if (ibin < 10000)
 633   //     bins[ibin] = bins [ibin-1]+1000;
 634   // }
 635
 636
 637   if(!h) {
 638     AliInfo(Form("Booking histo %s",hname.Data()));
 639     Bool_t oldStatus = TH1::AddDirectoryStatus();
 640     TH1::AddDirectory(kFALSE);
 641     //    h = new TH1F (hname.Data(), title, nbin, bins);
 642     //h = new TH1F (hname.Data(), title, 100, -0.5, 100-0.5);
 643     h = new TH1F (hname.Data(), title, 10000, -0.5, 10000-0.5);
 644     h->Sumw2();
 645     h->SetXTitle("SPD1");
 646     fHistoList->GetList()->Add(h);
 647     TH1::AddDirectory(oldStatus);
 648   }
 649   return h;
 650 }
 651
 652 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoV0M(const char * name, const char * title){
 653   // Book histo of events vs ntracklets, if needed
 654
 655   TString hname = "hV0M_";
 656   hname+=name;
 657   hname+=fHistoSuffix;
 658   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 659
 660   if(!h) {
 661     AliInfo(Form("Booking histo %s",hname.Data()));
 662     Bool_t oldStatus = TH1::AddDirectoryStatus();
 663     TH1::AddDirectory(kFALSE);
 664     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 665     h = new TH1F (hname.Data(), title, 500, -0.5, 500-0.5);
 666     h->Sumw2();
 667     h->SetXTitle("V0M");
 668     fHistoList->GetList()->Add(h);
 669     TH1::AddDirectory(oldStatus);
 670   }
 671   return h;
 672 }
 673
 674
 675 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracks(const char * name, const char * title){
 676   // Book histo of events vs ntracklets, if needed
 677
 678   TString hname = "hTPCTracks_";
 679   hname+=name;
 680   hname+=fHistoSuffix;
 681   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 682
 683   if(!h) {
 684     AliInfo(Form("Booking histo %s",hname.Data()));
 685     Bool_t oldStatus = TH1::AddDirectoryStatus();
 686     TH1::AddDirectory(kFALSE);
 687     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 688     h = new TH1F (hname.Data(), title, 50, -0.5, 50-0.5);
 689     h->Sumw2();
 690     h->SetXTitle("ntracks");
 691     fHistoList->GetList()->Add(h);
 692     TH1::AddDirectory(oldStatus);
 693   }
 694   return h;
 695 }
 696
 697
 698 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoEta(const char * name, const char * title){
 699   // Book histo of events vs ntracklets, if needed
 700
 701   TString hname = "hEta_";
 702   hname+=name;
 703   hname+=fHistoSuffix;
 704   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 705
 706   if(!h) {
 707     AliInfo(Form("Booking histo %s",hname.Data()));
 708     Bool_t oldStatus = TH1::AddDirectoryStatus();
 709     TH1::AddDirectory(kFALSE);
 710     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 711     h = new TH1F (hname.Data(), title, 20, -2, 2);
 712     h->Sumw2();
 713     h->SetXTitle("#eta");
 714     fHistoList->GetList()->Add(h);
 715     TH1::AddDirectory(oldStatus);
 716   }
 717   return h;
 718 }
 719
 720 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoPt(const char * name, const char * title){
 721   // Book histo of pt distribution, if needed
 722
 723   TString hname = "hPt_";
 724   hname+=name;
 725   hname+=fHistoSuffix;
 726   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 727
 728   if(!h) {
 729     AliInfo(Form("Booking histo %s",hname.Data()));
 730     Bool_t oldStatus = TH1::AddDirectoryStatus();
 731     TH1::AddDirectory(kFALSE);
 732     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 733     h = new TH1F (hname.Data(), title, 100, -0.5, 1-0.5);
 734     h->Sumw2();
 735     h->SetXTitle("p_{T}");
 736     fHistoList->GetList()->Add(h);
 737     TH1::AddDirectory(oldStatus);
 738   }
 739   return h;
 740 }
 741
 742 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationSPDTPCVz(const char * name, const char * title){
 743   // Book histo of events vs ntracklets, if needed
 744
 745   TString hname = "hTPCvsSPD_";
 746   hname+=name;
 747   hname+=fHistoSuffix;
 748   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 749
 750   if(!h) {
 751     AliInfo(Form("Booking histo %s",hname.Data()));
 752     Bool_t oldStatus = TH1::AddDirectoryStatus();
 753     TH1::AddDirectory(kFALSE);
 754     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 755     h = new TH2F (hname.Data(), title, 80, -20, 20,  80, -20, 20);
 756     h->Sumw2();
 757     h->SetXTitle("SPD Vz");
 758     h->SetYTitle("TPC Vz");
 759     fHistoList->GetList()->Add(h);
 760     TH1::AddDirectory(oldStatus);
 761   }
 762   return h;
 763 }
 764
 765 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationContrTPCSPDCls(const char * name, const char * title){
 766   // Book histo of events vs ntracklets, if needed
 767
 768   TString hname = "hContrTPCvsSPDCls_";
 769   hname+=name;
 770   hname+=fHistoSuffix;
 771   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 772
 773   if(!h) {
 774     AliInfo(Form("Booking histo %s",hname.Data()));
 775     Bool_t oldStatus = TH1::AddDirectoryStatus();
 776     TH1::AddDirectory(kFALSE);
 777     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 778     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
 779     h->Sumw2();
 780     h->SetXTitle("SPD Clusters");
 781     h->SetYTitle("TPC V Contributors");
 782     fHistoList->GetList()->Add(h);
 783     TH1::AddDirectory(oldStatus);
 784   }
 785   return h;
 786 }
 787 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationTrackletsSPDCls(const char * name, const char * title){
 788   // Book histo of events vs ntracklets, if needed
 789
 790   TString hname = "hTrackletsvsSPDCls_";
 791   hname+=name;
 792   hname+=fHistoSuffix;
 793   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 794
 795   if(!h) {
 796     AliInfo(Form("Booking histo %s",hname.Data()));
 797     Bool_t oldStatus = TH1::AddDirectoryStatus();
 798     TH1::AddDirectory(kFALSE);
 799     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 800     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
 801     h->Sumw2();
 802     h->SetXTitle("SPD Clusters");
 803     h->SetYTitle("N tracklets");
 804     fHistoList->GetList()->Add(h);
 805     TH1::AddDirectory(oldStatus);
 806   }
 807   return h;
 808 }
 809
 810
 811 void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const char * title, Bool_t * vdArray){
 812   //Fills a histo with the different trigger statistics in a venn like diagramm. Books it if needed.
 813
 814   // Get or book histo
 815   TString hname = "hTrigStat_";
 816   hname+=name;
 817   hname+=fHistoSuffix;
 818   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 819
 820   if(!h) {
 821     AliInfo(Form("Booking histo %s",hname.Data()));
 822     Bool_t oldStatus = TH1::AddDirectoryStatus();
 823     TH1::AddDirectory(kFALSE);
 824     Int_t nbins = 0;
 825     for(Int_t ientry = 0; ientry < kNVDEntries; ientry++){
 826       nbins = nbins | (1<<ientry);
 827     }
 828     //    cout << "NBINS " << nbins << endl;
 829     nbins = nbins+1;
 830     h = new TH1I (hname, title, nbins, -0.5, nbins-0.5);
 831     fHistoList->GetList()->Add(h);
 832     TH1::AddDirectory(oldStatus);
 833
 834     // we look at the combinations of n triggers
 835     // We set a bit for each trigger to fill the diagram
 836     // This is much simpler and faster than any recursive function
 837     h->GetXaxis()->SetBinLabel(1,"NONE");
 838     for(Int_t ibin = 1; ibin < nbins; ibin++){
 839       TString binname = "";
 840       Bool_t first = kTRUE;
 841       for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
 842         if (ibin & (1<<ivdentry)) {
 843           if(!first) binname += " & ";
 844           binname += kVDNames[ivdentry];
 845           first=kFALSE;
 846         }
 847       }
 848       h->GetXaxis()->SetBinLabel(ibin+1,binname.Data());
 849     }
 850
 851   }
 852
 853   UInt_t mask = 0;
 854   for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
 855     if(vdArray[ivdentry]) {
 856       mask  = mask | (1<<ivdentry);
 857       //      cout << " 1 "   ;
 858     } //else cout << " 0 ";
 859   }
 860   //  cout << hex << " = " << mask << endl;
 861
 862   h->Fill(mask);
 863
 864 }