1 /*************************************************************************
2 * Copyright(c) 1998-2008,ALICE Experiment at CERN,All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use,copy,modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee,provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /////////////////////////////////////////////////////
17 // AliAnalysisTaskFlowStrange:
18 // Analysis task to select K0/Lambda candidates for flow analysis.
19 // Authors: Cristian Ivan (civan@cern.ch)
20 // Carlos Perez (cperez@cern.ch)
21 // Pawel Debski (pdebski@cern.ch)
22 //////////////////////////////////////////////////////
31 #include "TProfile2D.h"
33 #include "TStopwatch.h"
38 #include "AliAnalysisManager.h"
39 #include "AliInputEventHandler.h"
41 #include "AliVVertex.h"
42 #include "AliVVZERO.h"
44 #include "AliMCEvent.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliESDVertex.h"
50 #include "AliESDtrackCuts.h"
52 #include "AliAODEvent.h"
53 #include "AliAODTrack.h"
54 #include "AliAODVertex.h"
56 #include "AliAODTracklets.h"
57 #include "AliAODHeader.h"
59 #include "AliAODMCHeader.h"
60 #include "AliAODMCParticle.h"
61 #include "TClonesArray.h"
62 #include "TDatabasePDG.h"
63 #include "TParticlePDG.h"
66 #include "TObjArray.h"
67 #include "AliFlowCandidateTrack.h"
69 #include "AliFlowTrackCuts.h"
70 #include "AliFlowEventCuts.h"
71 #include "AliFlowEvent.h"
72 #include "AliFlowBayesianPID.h"
73 #include "AliFlowCommonConstants.h"
74 #include "AliFlowVector.h"
76 #include "AliAnalysisTaskFlowStrange.h"
78 ClassImp(AliAnalysisTaskFlowStrange)
80 //=======================================================================
81 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange() :
95 fAddPiToMCReactionPlane(kTRUE),
98 fSkipSelection(kFALSE),
103 fExtraEventRejection(kFALSE),
104 fSkipCentralitySelection(kFALSE),
105 fCentMethod("V0MTRK"),
118 fExcludeTPCEdges(kFALSE),
153 fSkipTerminate(kTRUE),
173 fDecayDCAdaughters(0.0),
174 fDecayCosinePointingAngleXY(0.0),
176 fDecayDecayLength(0.0),
177 fDecayDecayLengthLab(0.0),
181 fDecayProductIPXY(0.0),
189 fDecayMatchOrigin(0.0),
193 fDecayMatchRadXY(0.0),
197 fDecayMaxDCAdaughters(0.0),
198 fDecayMinCosinePointingAngleXY(0.0),
200 fDecayAPCutPie(kTRUE),
201 fDecayStopPIDAtPt(3.0),
203 fDecayMaxDecayLength(0.0),
204 fDecayMaxProductIPXY(0.0),
205 fDecayMaxRapidity(0.0),
212 fDaughterNFClsTPC(0),
213 fDaughterNSClsTPC(0),
214 fDaughterChi2PerNClsTPC(0.0),
216 fDaughterImpactParameterXY(0.0),
217 fDaughterImpactParameterZ(0.0),
220 fDaughterNSigmaPID(0.0),
221 fDaughterKinkIndex(0),
222 fDaughterAtSecPhi(0.0),
223 fDaughterAtSecEta(0.0),
224 fDaughterAtSecPt(0.0),
225 fDaughterMatchPhi(0.0),
226 fDaughterMatchEta(0.0),
227 fDaughterMatchPt(0.0),
228 fDaughterMatchImpactParameterXY(0.0),
229 fDaughterMatchImpactParameterZ(0.0),
230 fDaughterUnTag(kTRUE),
231 fDaughterMinEta(0.0),
232 fDaughterMaxEta(0.0),
234 fDaughterMinNClsTPC(0),
235 fDaughterMinNClsITS(-1),
236 fDaughterMinXRows(0),
237 fDaughterMaxChi2PerNClsTPC(0.0),
238 fDaughterMinXRowsOverNClsFTPC(0.0),
239 fDaughterMinImpactParameterXY(0.0),
240 fDaughterMaxNSigmaPID(0.0),
241 fDaughterSPDRequireAny(kFALSE),
242 fDaughterITSrefit(kFALSE) {
244 for(Int_t i=0; i!=100; ++i) fPtBinEdge[i]=0;
245 for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
246 for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
247 for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
248 for(Int_t i=0; i!=64; ++i) fVZEextW[i]=1;
250 //=======================================================================
251 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange(const char *name) :
252 AliAnalysisTaskSE(name),
265 fAddPiToMCReactionPlane(kTRUE),
268 fSkipSelection(kFALSE),
273 fExtraEventRejection(kFALSE),
274 fSkipCentralitySelection(kFALSE),
275 fCentMethod("V0MTRK"),
288 fExcludeTPCEdges(kFALSE),
323 fSkipTerminate(kTRUE),
343 fDecayDCAdaughters(0.0),
344 fDecayCosinePointingAngleXY(0.0),
346 fDecayDecayLength(0.0),
347 fDecayDecayLengthLab(0.0),
351 fDecayProductIPXY(0.0),
359 fDecayMatchOrigin(0.0),
363 fDecayMatchRadXY(0.0),
367 fDecayMaxDCAdaughters(0.0),
368 fDecayMinCosinePointingAngleXY(0.0),
370 fDecayAPCutPie(kTRUE),
371 fDecayStopPIDAtPt(3.0),
373 fDecayMaxDecayLength(0.0),
374 fDecayMaxProductIPXY(0.0),
375 fDecayMaxRapidity(0.0),
382 fDaughterNFClsTPC(0),
383 fDaughterNSClsTPC(0),
384 fDaughterChi2PerNClsTPC(0.0),
386 fDaughterImpactParameterXY(0.0),
387 fDaughterImpactParameterZ(0.0),
390 fDaughterNSigmaPID(0.0),
391 fDaughterKinkIndex(0),
392 fDaughterAtSecPhi(0.0),
393 fDaughterAtSecEta(0.0),
394 fDaughterAtSecPt(0.0),
395 fDaughterMatchPhi(0.0),
396 fDaughterMatchEta(0.0),
397 fDaughterMatchPt(0.0),
398 fDaughterMatchImpactParameterXY(0.0),
399 fDaughterMatchImpactParameterZ(0.0),
400 fDaughterUnTag(kTRUE),
401 fDaughterMinEta(0.0),
402 fDaughterMaxEta(0.0),
404 fDaughterMinNClsTPC(0),
405 fDaughterMinNClsITS(-1),
406 fDaughterMinXRows(0),
407 fDaughterMaxChi2PerNClsTPC(0.0),
408 fDaughterMinXRowsOverNClsFTPC(0.0),
409 fDaughterMinImpactParameterXY(0.0),
410 fDaughterMaxNSigmaPID(0.0),
411 fDaughterSPDRequireAny(kFALSE),
412 fDaughterITSrefit(kFALSE) {
414 for(Int_t i=0; i!=100; ++i) fPtBinEdge[i]=0;
415 for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
416 for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
417 for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
418 for(Int_t i=0; i!=64; ++i) fVZEextW[i]=1;
419 DefineInput( 0,TChain::Class());
420 DefineOutput(1,TList::Class());
421 DefineOutput(2,AliFlowEventSimple::Class()); // TPC object
422 DefineOutput(3,AliFlowEventSimple::Class()); // VZE object
424 //=======================================================================
425 AliAnalysisTaskFlowStrange::~AliAnalysisTaskFlowStrange() {
427 if (fCandidates) delete fCandidates;
428 if (fTPCevent) delete fTPCevent;
429 if (fVZEevent) delete fVZEevent;
430 if (fList) delete fList;
432 //=======================================================================
433 void AliAnalysisTaskFlowStrange::SetPtEdges(Int_t n, Double_t *p) {
435 for(int i=0;i!=n+1;++i) fPtBinEdge[i] = p[i];
437 //=======================================================================
438 TList* AliAnalysisTaskFlowStrange::RunTerminateAgain(TList *lst) {
439 if(!lst) return NULL;
441 fSpecie = Int_t( ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSpecie) );
442 fSkipSelection = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSkipSelection);
443 fReadMC = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kReadMC);
447 //=======================================================================
448 void AliAnalysisTaskFlowStrange::PrintConfig() {
450 printf("******************************\n");
451 printf("<TASK Configuration> %s\n",GetName());
452 printf(" fDebug %d\n",fDebug);
453 printf(" fQAlevel %d\n",fQAlevel);
454 printf(" fExtraEventRejection %s\n",fExtraEventRejection?"kTRUE":"kFALSE");
455 printf(" fCentMethod %s\n",fCentMethod.Data());
456 printf(" fCentPerMin %d\n",fCentPerMin);
457 printf(" fCentPerMax %d\n",fCentPerMax);
458 printf(" fVextexZcut %f\n",fVertexZcut);
459 printf(" fRunOnpA %s\n",fRunOnpA?"kTRUE":"kFALSE");
460 printf(" fRunOnpp %s\n",fRunOnpp?"kTRUE":"kFALSE");
461 printf(" fReadESD %s\n",fReadESD?"kTRUE":"kFALSE");
462 printf(" fReadMC %s\n",fReadMC?"kTRUE":"kFALSE");
464 printf(" fAddPiToMCReactionPlane %s\n",fAddPiToMCReactionPlane?"kTRUE":"kFALSE");
465 printf(" fPostMatched %d\n",fPostMatched);
466 printf(" fAvoidExec %s\n",fAvoidExec?"kTRUE":"kFALSE");
467 printf(" fSkipCentralitySelection %s\n",fSkipCentralitySelection?"kTRUE":"kFALSE");
469 printf(" fVZEsave %s\n",fVZEsave?"kTRUE":"kFALSE");
471 printf(" fVZEload %d runs\n",fVZEload->GetEntries());
472 printf(" fVZEmb %s\n",fVZEmb?"kTRUE":"kFALSE");
473 printf(" fVZEByDisk %s\n",fVZEByDisk?"kTRUE":"kFALSE");
475 printf(" fHarmonic %d\n",fHarmonic);
476 printf(" fWhichPsi %d\n",fWhichPsi);
477 printf(" fVZECa %d\n",fVZECa);
478 printf(" fVZECb %d\n",fVZECb);
479 printf(" fVZEAa %d\n",fVZEAa);
480 printf(" fVZEAb %d\n",fVZEAb);
481 printf(" fRFPFilterBit %d\n",fRFPFilterBit);
482 printf(" fRFPminPt %f\n",fRFPminPt);
483 printf(" fRFPmaxPt %f\n",fRFPmaxPt);
484 printf(" fRFPAminEta %f\n",fRFPAminEta);
485 printf(" fRFPAmaxEta %f\n",fRFPAmaxEta);
486 printf(" fRFPCminEta %f\n",fRFPCminEta);
487 printf(" fRFPCmaxEta %f\n",fRFPCmaxEta);
488 printf(" fRFPmaxIPxy %f\n",fRFPmaxIPxy);
489 printf(" fRFPmaxIPz %f\n",fRFPmaxIPz);
490 printf(" fRFPTPCsignal %f\n",fRFPTPCsignal);
491 printf(" fRFPTPCncls %d\n",fRFPTPCncls);
492 printf(" fExcludeTPCEdges %s\n",fExcludeTPCEdges?"kTRUE":"kFALSE");
493 printf(" fSkipSelection %s\n",fSkipSelection?"kTRUE":"kFALSE");
494 if(!fSkipSelection) {
495 printf(" fSpecie %d\n",fSpecie);
496 printf(" fPtBins %d\n |",fPtBins);
497 for(int i=0; i!=fPtBins+1; ++i) printf("%f|",fPtBinEdge[i]); printf("\n");
499 printf(" fMassBins %d\n",fMassBins);
500 printf(" fMinMass %f\n",fMinMass);
501 printf(" fMaxMass %f\n",fMaxMass);
504 printf(" fSkipVn %s\n",fSkipVn?"kTRUE":"kFALSE");
506 printf(" fUseFP %s\n",fUseFP?"kTRUE":"kFALSE");
510 //=======================================================================
511 void AliAnalysisTaskFlowStrange::MyPrintConfig() {
512 // Dump for derived task
513 printf("==================================\n");
514 printf("<FlowStrange> \n");
515 if(!fSkipSelection) {
517 printf(" fOnline %s\n",fOnline?"kTRUE":"kFALSE");
518 printf(" fHomemade %s\n",fHomemade?"kTRUE":"kFALSE");
520 printf(" fDecayMinEta %f\n",fDecayMinEta);
521 printf(" fDecayMaxEta %f\n",fDecayMaxEta);
522 printf(" fDecayMinPt %f\n",fDecayMinPt);
523 printf(" fDecayMaxDCAdaughters %f\n",fDecayMaxDCAdaughters);
524 printf(" fDecayMinCosinePointingAngleXY %f\n",fDecayMinCosinePointingAngleXY);
525 printf(" fDecayMinQt %f\n",fDecayMinQt);
526 printf(" fDecayAPCutPie %s\n",fDecayAPCutPie?"kTRUE":"kFALSE");
527 printf(" fDecayStopPIDAtPt %f\n",fDecayStopPIDAtPt);
528 printf(" fDecayMinRadXY %f\n",fDecayMinRadXY);
529 printf(" fDecayMaxDecayLength %f\n",fDecayMaxDecayLength);
530 printf(" fDecayMaxProductIPXY %f\n",fDecayMaxProductIPXY);
531 printf(" fDecayMaxRapidity %f\n",fDecayMaxRapidity);
533 printf(" fDaughterUnTag %s\n",fDaughterUnTag?"kTRUE":"kFALSE");
534 printf(" fDaughterMinEta %f\n",fDaughterMinEta);
535 printf(" fDaughterMaxEta %f\n",fDaughterMaxEta);
536 printf(" fDaughterMinPt %f\n",fDaughterMinPt);
537 printf(" fDaughterMinNClsTPC %d\n",fDaughterMinNClsTPC);
538 printf(" fDaughterMinXRows %d\n",fDaughterMinXRows);
539 printf(" fDaughterMaxChi2PerNClsTPC %f\n",fDaughterMaxChi2PerNClsTPC);
540 printf(" fDaughterMinXRowsOverNClsFTPC %f\n",fDaughterMinXRowsOverNClsFTPC);
541 printf(" fDaughterMinImpactParameterXY %f\n",fDaughterMinImpactParameterXY);
542 printf(" fDaughterMaxNSigmaPID %f\n",fDaughterMaxNSigmaPID);
544 //=======================================================================
545 void AliAnalysisTaskFlowStrange::UserCreateOutputObjects() {
546 //UserCreateOutputObjects
547 if(fDebug) PrintConfig();
552 if(fReadESD) MakeFilterBits();
554 AliFlowCommonConstants *cc = AliFlowCommonConstants::GetMaster();
555 cc->SetNbinsMult(3000); cc->SetMultMin(0); cc->SetMultMax(30000);
556 cc->SetNbinsPt(100); cc->SetPtMin(0.0); cc->SetPtMax(20.0);
557 cc->SetNbinsPhi(100); cc->SetPhiMin(0.0); cc->SetPhiMax(TMath::TwoPi());
558 cc->SetNbinsEta(100); cc->SetEtaMin(-5.0); cc->SetEtaMax(+5.0);
559 cc->SetNbinsQ(100); cc->SetQMin(0.0); cc->SetQMax(3.0);
560 cc->SetNbinsMass(fMassBins);
561 cc->SetMassMin(fMinMass);
562 cc->SetMassMax(fMaxMass);
564 //loading pid response
565 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
566 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
567 fPIDResponse = inputHandler->GetPIDResponse();
570 fTPCevent = new AliFlowEvent(100);
571 fVZEevent = new AliFlowEvent(100);
572 //array of candidates
573 fCandidates = new TObjArray(100);
574 fCandidates->SetOwner();
577 if(fUseFP) { // for connection to the flow package
578 PostData(2,fTPCevent);
579 PostData(3,fVZEevent);
584 //=======================================================================
585 void AliAnalysisTaskFlowStrange::MyUserCreateOutputObjects() {
592 tList=new TList(); tList->SetName("ESD_TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
593 tList=new TList(); tList->SetName("ESD_TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
594 tH2D = new TH2D("NPAIR", "NPAIR;NPOS;NNEG",1000,0,5000,1000,0,5000); tList->Add(tH2D);
595 tH2D = new TH2D("PtIPXY","PtIPXY;Pt;IPxy", 100,0,10,200,-10,+10); tList->Add(tH2D);
597 //aod prefilter candidates
598 tList=new TList(); tList->SetName("V0SAll"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
599 tH2D = new TH2D("V0SADC","V0S AFTER DAUGHTER CUTS;V0ALL;V0IMW",100,0,1000,100,0,1000); tList->Add(tH2D);
600 tList=new TList(); tList->SetName("AllDau"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
602 tList=new TList(); tList->SetName("V0SSel"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
603 tList=new TList(); tList->SetName("SelDau"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
606 tList=new TList(); tList->SetName("V0SAllVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
607 tList=new TList(); tList->SetName("V0SSelVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
611 tList=new TList(); tList->SetName("V0SAllIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
612 tList=new TList(); tList->SetName("V0SAllOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
613 tList=new TList(); tList->SetName("V0SSelIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
614 tList=new TList(); tList->SetName("V0SSelOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
618 tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
619 tH1D = new TH1D("Events", "Events",5,0.5,5.5); tList->Add(tH1D);
620 tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
621 tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
622 tH1D->GetXaxis()->SetBinLabel(3,"Daughters in stack");
623 tH1D->GetXaxis()->SetBinLabel(4,"Correspond to decay");
624 tH1D->GetXaxis()->SetBinLabel(5,"Decay has mother");
625 tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
626 tList=new TList(); tList->SetName("MthDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
627 tList=new TList(); tList->SetName("MthPosPos"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
628 tList=new TList(); tList->SetName("MthNegNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
629 tList=new TList(); tList->SetName("MthPosNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
630 tList=new TList(); tList->SetName("MthNegDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
631 tList=new TList(); tList->SetName("MthPosDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
632 tList=new TList(); tList->SetName("MthFeedDown"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
633 tList=new TList(); tList->SetName("UnMth"); tList->SetOwner(); AddCandidatesSpy(tList,false); fList->Add(tList);
634 tList=new TList(); tList->SetName("UnMthDau"); tList->SetOwner(); AddTrackSpy(tList,false); fList->Add(tList);
636 tList=new TList(); tList->SetName("V0SMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
637 tList=new TList(); tList->SetName("V0SMthPosPosVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
638 tList=new TList(); tList->SetName("V0SMthNegNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
639 tList=new TList(); tList->SetName("V0SMthPosNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
640 tList=new TList(); tList->SetName("V0SUnMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
644 //=======================================================================
645 void AliAnalysisTaskFlowStrange::AddQAEvents() {
646 // function to add event qa
649 TList *tQAEvents=new TList();
650 tQAEvents->SetName("Event");
651 tQAEvents->SetOwner();
652 tH1D = new TH1D("Events","Number of Events",6,0,6); tQAEvents->Add(tH1D);
653 tH1D->GetXaxis()->SetBinLabel(1,"exec");
654 tH1D->GetXaxis()->SetBinLabel(2,"userexec");
655 tH1D->GetXaxis()->SetBinLabel(3,"reached");
656 tH1D->GetXaxis()->SetBinLabel(4,"selected");
657 tH1D->GetXaxis()->SetBinLabel(5,"rejectedByLowQw");
658 tH1D->GetXaxis()->SetBinLabel(6,"rejectedByErrorLoadVZEcal");
659 tProfile = new TProfile("Configuration","Configuration",10,0.5,10.5); tQAEvents->Add(tProfile);
660 tProfile->Fill(kSpecie,fSpecie,1);
661 tProfile->GetXaxis()->SetBinLabel(kSpecie,"fSpecie");
662 tProfile->Fill(kHarmonic,fHarmonic,1);
663 tProfile->GetXaxis()->SetBinLabel(kHarmonic,"fHarmonic");
664 tProfile->Fill(kReadMC,fReadMC,1);
665 tProfile->GetXaxis()->SetBinLabel(kReadMC,"fReadMC");
666 tProfile->Fill(kSkipSelection,fSkipSelection,1);
667 tProfile->GetXaxis()->SetBinLabel(kSkipSelection,"fSkipSelection");
668 tH1D = new TH1D("POI","POIs;multiplicity",800,0,800); tQAEvents->Add(tH1D);
669 tH1D = new TH1D("UNTAG","UNTAG;Untagged Daughters",800,0,800);tQAEvents->Add(tH1D);
670 tH1D = new TH1D("RealTime","RealTime;LogT sec",2000,-10,+10); tQAEvents->Add(tH1D);
671 fList->Add(tQAEvents);
672 AddEventSpy("EventsRaw");
673 AddEventSpy("EventsReached");
674 AddEventSpy("EventsSelected");
675 AddEventSpy("EventsAnalyzed");
679 //=======================================================================
680 void AliAnalysisTaskFlowStrange::AddEventSpy(TString name) {
683 TList *tList=new TList();
684 tList->SetName(name.Data());
686 tH2D = new TH2D("VTXZ","VTXZ;PriVtxZ;SPDVtxZ",60,-25,+25,60,-25,+25); tList->Add( tH2D );
687 tH2D = new TH2D("CCCC","CCCC;V0M;TRK",60,-10,110,60,-10,110); tList->Add( tH2D );
688 tH2D = new TH2D("HYBTPC","HYBTPC;TPC ONLY;HYBRID",100,0,3000,100,0,3000); tList->Add( tH2D );
689 tH1D = new TH1D("HYBTPCRat","HYBTPCRat;TPC/HYB",120,0.2,2.2); tList->Add( tH1D );
690 tH2D = new TH2D("SPDVZE","SPDVZE;SPD Tracklets;Total Multiplicity in VZERO",100,0,3500,100,0,25000); tList->Add( tH2D );
691 tH1D = new TH1D("SPDVZERat","SPDVZERat;TotalMultiplicityVZERO/SPDTracklets",120,2,+12); tList->Add( tH1D );
693 tH1D = new TH1D("MCEP","MCEP;MCEP",100,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
697 //=======================================================================
698 void AliAnalysisTaskFlowStrange::FillEventSpy(TString name) {
699 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("VTXZ"))->Fill( fPriVtxZ, fSPDVtxZ );
700 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("CCCC"))->Fill( fV0M, fTRK );
701 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("HYBTPC"))->Fill( fRefMultTPC, fRefMultHyb );
703 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("HYBTPCRat"))->Fill( double(fRefMultTPC)/double(fRefMultHyb) );
704 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("SPDVZE"))->Fill( fSPDtracklets, fVZETotM );
706 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("SPDVZERat"))->Fill( fVZETotM/fSPDtracklets );
708 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("MCEP"))->Fill( fMCEP );
711 //=======================================================================
712 void AliAnalysisTaskFlowStrange::AddMakeQSpy() {
716 TList *tList=new TList();
717 tList->SetName("MakeQSpy");
720 tH1D = new TH1D("RFPTPC","TPC Refrence Multiplicity;multiplicity",3000,0,3000); tList->Add( tH1D );
721 tH1D = new TH1D("RFPVZE","VZERO Reference Multiplicity;multiplicity",3000,0,30000); tList->Add( tH1D );
722 tH1D = new TH1D("QmTPC","TPC Normalized Q vector;|Q|/#sqrt{M}",360,0,7); tList->Add( tH1D );
723 tH1D = new TH1D("QmVZEA","VZEROA Normalized Q vector;|Q|/#sqrt{W}",360,0,7); tList->Add( tH1D );
724 tH1D = new TH1D("QmVZEC","VZEROC Normalized Q vector;|Q|/#sqrt{W}",360,0,7); tList->Add( tH1D );
725 tH2D = new TH2D("TPCAllPhiEta","TPCall;Phi;Eta",180,0,TMath::TwoPi(),80,-0.9,+0.9); tList->Add( tH2D );
726 tH2D = new TH2D("VZEAllPhiEta","VZEall;Phi;Eta",20,0,TMath::TwoPi(),40,-4.0,+6.0); tList->Add( tH2D );
727 tH1D = new TH1D("TPCPSI","TPCPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
728 tH1D = new TH1D("TPCPSIA","TPCPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
729 tH1D = new TH1D("TPCPSIC","TPCPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
730 tH1D = new TH1D("VZEPSI","VZEPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
731 tH1D = new TH1D("VZEPSIA","VZEPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
732 tH1D = new TH1D("VZEPSIC","VZEPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
733 tH2D = new TH2D("PSI_TPCAVZEC","PSI_TPCAVZEC",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
734 tH2D = new TH2D("PSI_TPCCVZEA","PSI_TPCAVZEC",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
735 tH2D = new TH2D("PSI_TPCVZE","PSI_TPCVZE",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
736 tPF1 = new TProfile("TPCQm","TPCQm",6,0.5,6.5); tList->Add( tPF1 );
737 tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
738 tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
739 tPF1->GetXaxis()->SetBinLabel(5,"Qy"); tPF1->GetXaxis()->SetBinLabel(6,"Qx");
740 tPF1 = new TProfile("VZEQm","VZEQm",6,0.5,6.5); tList->Add( tPF1 );
741 tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
742 tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
743 tPF1->GetXaxis()->SetBinLabel(5,"Qy"); tPF1->GetXaxis()->SetBinLabel(6,"Qx");
744 tPF1 = new TProfile("QmVZEAQmVZEC","QmVZEAQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
745 tPF1 = new TProfile("QmVZEASQUARED","QmVZEASQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
746 tPF1 = new TProfile("QmVZECSQUARED","QmVZECSQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
747 tPF1 = new TProfile("QmTPCQmVZEA","QmTPCQmVZEA",1,0.5,1.5,"s"); tList->Add( tPF1 );
748 tPF1 = new TProfile("QmTPCQmVZEC","QmTPCQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
749 tH1D = new TH1D("ChiSquaredVZEA","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
750 tH1D = new TH1D("ChiSquaredVZEC","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
752 tH1D = new TH1D("PSIMCDIFFTPC","PSIMCDIFFTPC;MC-TPC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
753 tH1D = new TH1D("PSIMCDIFFTPCA","PSIMCDIFFTPCA;MC-TPCA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
754 tH1D = new TH1D("PSIMCDIFFTPCC","PSIMCDIFFTPCC;MC-TPCC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
755 tH1D = new TH1D("PSIMCDIFFVZE","PSIMCDIFFVZE;MC-VZE",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
756 tH1D = new TH1D("PSIMCDIFFVZEA","PSIMCDIFFVZEA;MC-VZEA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
757 tH1D = new TH1D("PSIMCDIFFVZEC","PSIMCDIFFVZEC;MC-VZEC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
759 tList=new TList(); tList->SetName("TPCRFPall"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
760 tList=new TList(); tList->SetName("TPCRFPsel"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
762 //=======================================================================
763 void AliAnalysisTaskFlowStrange::AddQACandidates() {
764 // function to add histogramming for candidates
765 if(fSkipSelection) return;
769 //charge particles (benchmark)
771 tList=new TList(); tList->SetName("TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
772 tList=new TList(); tList->SetName("TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
774 tList=new TList(); tList->SetName("TrkAllVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
775 tList=new TList(); tList->SetName("TrkSelVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
779 tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
780 tH1D = new TH1D("Events", "Events",3,0.5,3.5); tList->Add(tH1D);
781 tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
782 tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
783 tH1D->GetXaxis()->SetBinLabel(3,"Track in stack");
784 tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
785 tList=new TList(); tList->SetName("MthPos"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
786 tList=new TList(); tList->SetName("MthNeg"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
788 tList=new TList(); tList->SetName("MthVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
789 tList=new TList(); tList->SetName("MthPosVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
790 tList=new TList(); tList->SetName("MthNegVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
796 tList=new TList(); tList->SetName("MCTPionGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
797 tList=new TList(); tList->SetName("MCTKaonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
798 tList=new TList(); tList->SetName("MCTK0sGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
799 tList=new TList(); tList->SetName("MCTProtonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
800 tList=new TList(); tList->SetName("MCTLdaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
801 tList=new TList(); tList->SetName("MCTPhiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
802 tList=new TList(); tList->SetName("MCTXiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
803 tList=new TList(); tList->SetName("MCTOmegaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
804 tList=new TList(); tList->SetName("MCTPion"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
805 tList=new TList(); tList->SetName("MCTKaon"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
806 tList=new TList(); tList->SetName("MCTK0s"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
807 tList=new TList(); tList->SetName("MCTLda"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
808 tList=new TList(); tList->SetName("MCTProton"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
810 MyUserCreateOutputObjects();
812 //=======================================================================
813 void AliAnalysisTaskFlowStrange::Exec(Option_t* option) {
814 // bypassing ::exec (needed because of AMPT)
815 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(0);
817 AliAnalysisTaskFlowStrange::UserExec(option);
819 AliAnalysisTaskSE::Exec(option);
822 //=======================================================================
823 void AliAnalysisTaskFlowStrange::UserExec(Option_t *option) {
825 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(1);
826 AliAnalysisTaskFlowStrange::MyUserExec(option);
828 //=======================================================================
829 void AliAnalysisTaskFlowStrange::MyNotifyRun() {
830 if(fVZEsave) AddVZEROResponse();
832 //=======================================================================
833 Bool_t AliAnalysisTaskFlowStrange::CalibrateEvent() {
834 if(fVZEsave) SaveVZEROResponse();
838 if(!fVZEResponse) okay = kFALSE;
842 //=======================================================================
843 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliESDEvent*) {
844 // ESD reading discontinued: TO BE UPDATED
846 Double_t acceptEvent=kTRUE;
847 Double_t tTPCVtxZ = tESD->GetPrimaryVertexTPC()->GetZ();
848 if(tESD->GetPrimaryVertexTPC()->GetNContributors()<=0) return kFALSE;
849 Double_t tSPDVtxZ = tESD->GetPrimaryVertexSPD()->GetZ();
850 if(tESD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
852 AliCentrality *cent = tESD->GetCentrality();
854 cc1 = cent->GetCentralityPercentile("V0M");
855 cc2 = cent->GetCentralityPercentile("TRK");
856 TString mycent = fCentMethod;
857 if(fCentMethod.Contains("V0MTRK")) {
858 acceptEvent = TMath::Abs(cc1-cc2)>5.0?kFALSE:acceptEvent; // a la Alex
861 fThisCent = cent->GetCentralityPercentile( mycent );
862 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
863 acceptEvent = TMath::Abs(tTPCVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
864 acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
865 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
866 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
869 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
870 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
876 //=======================================================================
877 Bool_t AliAnalysisTaskFlowStrange::MinimumRequirementsAA(AliAODEvent *tAOD) {
878 fRunNumber = tAOD->GetRunNumber();
879 AliCentrality *cent = ((AliVAODHeader*)tAOD->GetHeader())->GetCentralityP();
880 fV0M = cent->GetCentralityPercentile("V0M");
881 fTRK = cent->GetCentralityPercentile("TRK");
882 TString mycent = fCentMethod;
883 if(fCentMethod.Contains("V0MTRK")) {
886 fThisCent = cent->GetCentralityPercentile( mycent );
887 fPriVtxZ = tAOD->GetPrimaryVertex()->GetZ();
888 fSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
889 fSPDtracklets = tAOD->GetTracklets()->GetNumberOfTracklets();
890 fVZETotM = tAOD->GetVZEROData()->GetMTotV0A() + tAOD->GetVZEROData()->GetMTotV0C();
891 int hyb_fb = 272; // for 2010h::AOD086
892 if(fRunNumber>=166529&&fRunNumber<=170593) {
893 hyb_fb = 768; // for 2011h::AOD145
895 fRefMultTPC = RefMult(tAOD,128);
896 fRefMultHyb = RefMult(tAOD,hyb_fb);
899 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(tAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
901 fMCEP = mcHeader->GetReactionPlaneAngle();
902 if(fAddPiToMCReactionPlane) fMCEP += (gRandom->Rndm()>0.5)*TMath::Pi();
905 // centrality selection health
906 // cut in Vtx 10 & NContributors
907 if(!fSkipCentralitySelection) if(fThisCent<0||fThisCent>100) return kFALSE;
908 // vtx z position compatibility within 5 mm
909 if(TMath::Abs(fPriVtxZ-fSPDVtxZ)>0.5) return kFALSE;
910 if(fExtraEventRejection) {
911 // specific cuts for 2010h (AOD086)
912 if(fRunNumber>=136851&&fRunNumber<=139517) {
913 if(fRefMultTPC>1.118*fRefMultHyb+100) return kFALSE;
914 if(fRefMultTPC<1.118*fRefMultHyb-100) return kFALSE;
916 // specific cuts for 2011h (AOD145)
917 if(fRunNumber>=166529&&fRunNumber<=170593) {
918 if(fRefMultTPC>1.205*fRefMultHyb+100) return kFALSE;
919 if(fRefMultTPC<1.205*fRefMultHyb-100) return kFALSE;
924 //=======================================================================
925 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliAODEvent *tAOD) {
926 Bool_t minimum = MinimumRequirementsAA(tAOD);
927 FillEventSpy("EventsRaw");
928 if(!minimum) return kFALSE;
930 Double_t acceptEvent=kTRUE;
931 TString mycent = fCentMethod;
932 if(fCentMethod.Contains("V0MTRK")) {
933 acceptEvent = TMath::Abs(fV0M-fTRK)>5.0?kFALSE:acceptEvent;
936 if(!fSkipCentralitySelection) acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
937 acceptEvent = TMath::Abs(fPriVtxZ)>fVertexZcut?kFALSE:acceptEvent;
939 FillEventSpy("EventsReached");
940 if(acceptEvent) FillEventSpy("EventsSelected");
943 //=======================================================================
944 Bool_t AliAnalysisTaskFlowStrange::AcceptPPEvent(AliAODEvent*) {
945 // PP reading discontinued: TO BE UPDATED
947 Double_t acceptEvent=kTRUE;
948 Double_t tVtxZ = tAOD->GetPrimaryVertex()->GetZ();
949 if(tAOD->GetPrimaryVertex()->GetNContributors()<=0) return kFALSE;
950 Double_t tSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
951 if(tAOD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
953 AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
955 cc1 = cent->GetCentralityPercentile("V0M");
956 cc2 = cent->GetCentralityPercentile("TRK");
957 fThisCent = GetReferenceMultiplicity();
958 //for pp i use fCentPerXXX to select on multiplicity
959 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
960 acceptEvent = TMath::Abs(tVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
961 acceptEvent = TMath::Abs(tVtxZ)>fVertexZcut?kFALSE:acceptEvent;
962 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
963 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
966 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
967 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
973 //=======================================================================
974 Int_t AliAnalysisTaskFlowStrange::GetReferenceMultiplicity() { //toberefined
975 AliAODEvent *tAOD = (AliAODEvent *) InputEvent();
978 Int_t rawN = tAOD->GetNumberOfTracks();
980 for(Int_t id=0; id!=rawN; ++id) {
981 track = dynamic_cast<AliAODTrack*>(tAOD->GetTrack(id));
982 if(!track) AliFatal("Not a standard AOD");
983 if(!track->TestFilterBit(fRFPFilterBit)) continue;
988 //=======================================================================
989 Bool_t AliAnalysisTaskFlowStrange::AcceptPAEvent(AliAODEvent*) {
990 // PA reading discontinued: TO BE UPDATED
992 //if(aod->GetHeader()->GetEventNumberESDFile() == 0) return; //rejecting first chunk NOT NEEDED ANYMORE
993 Int_t bc2 = ((AliVAODHeader*)tAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
994 if(bc2!=0) return kFALSE;
995 Int_t bc1 = ((AliVAODHeader*)tAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
996 if(bc1!=0) return kFALSE;
997 Short_t isPileup = tAOD->IsPileupFromSPD(5);
998 if(isPileup!=0) return kFALSE;
999 if(tAOD->GetHeader()->GetRefMultiplicityComb08()<0) return kFALSE;
1001 const AliAODVertex* spdVtx = tAOD->GetPrimaryVertexSPD();
1002 if(!spdVtx) return kFALSE;
1003 if(spdVtx->GetNContributors()<=0) return kFALSE;
1005 const AliAODVertex* tpcVtx=NULL;
1006 Int_t nVertices = tAOD->GetNumberOfVertices();
1007 for(Int_t iVertices = 0; iVertices < nVertices; iVertices++){
1008 const AliAODVertex* vertex = tAOD->GetVertex(iVertices);
1009 if (vertex->GetType() != AliAODVertex::kMainTPC) continue;
1012 if(!tpcVtx) return kFALSE;
1013 if(tpcVtx->GetNContributors()<=0) return kFALSE;
1014 Double_t tTPCVtxZ = tpcVtx->GetZ();
1015 Double_t tSPDVtxZ = spdVtx->GetZ();
1016 if (TMath::Abs(tSPDVtxZ - tTPCVtxZ)>2.0) return kFALSE;
1017 if(plpMV(tAOD)) return kFALSE;
1019 Double_t acceptEvent=kTRUE;
1021 AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
1023 cc1 = cent->GetCentralityPercentile("V0M");
1024 cc2 = cent->GetCentralityPercentile("TRK");
1025 if(fCentMethod.Contains("V0MTRK")) fCentMethod = "V0M";
1026 fThisCent = cent->GetCentralityPercentile( fCentMethod );
1027 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
1028 acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
1030 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
1031 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
1033 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
1034 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
1040 //=======================================================================
1041 void AliAnalysisTaskFlowStrange::MyUserExec(Option_t *) {
1046 printf("****************\n");
1047 printf("****************\n");
1048 printf("**::MyUserExec()\n");
1050 if(fUseFP) fCandidates->SetLast(-1);
1051 AliESDEvent *tESD=dynamic_cast<AliESDEvent*>(InputEvent());
1052 AliAODEvent *tAOD=dynamic_cast<AliAODEvent*>(InputEvent());
1053 Int_t prevRun = fRunNumber;
1055 Bool_t acceptEvent=kFALSE;
1057 if(!tESD) {ResetContainers(); Publish(); return;}
1058 acceptEvent = fRunOnpp?kFALSE:fRunOnpA?kFALSE:AcceptAAEvent(tESD);
1060 if(!tAOD) {ResetContainers(); Publish(); return;}
1061 acceptEvent = fRunOnpp?AcceptPPEvent(tAOD):fRunOnpA?AcceptPAEvent(tAOD):AcceptAAEvent(tAOD);
1063 if(prevRun!=fRunNumber) {
1066 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(2);
1067 //=>does the event clear?
1068 if(!acceptEvent) {ResetContainers(); Publish(); return;}
1069 // healthy event incomming
1070 if( !CalibrateEvent() ) { // saves/retrieves/qas VZEROCAL
1071 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(5);
1072 ResetContainers(); Publish(); return; // issue retrieving callibration
1077 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(4);
1078 ResetContainers(); Publish(); return;
1081 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(3);
1082 //=>great, lets do our stuff!
1083 FillEventSpy("EventsAnalyzed");
1086 if(!fSkipSelection) {
1090 if(fSpecie<10) ReadFromAODv0(tAOD);
1091 else ChargeParticles(tAOD);
1093 if(fUseFP) AddCandidates();
1098 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("RealTime"))->Fill( TMath::Log( tTime.RealTime() ) );
1101 //=======================================================================
1102 void AliAnalysisTaskFlowStrange::Publish() {
1105 PostData(2,fTPCevent);
1106 PostData(3,fVZEevent);
1109 //=======================================================================
1110 void AliAnalysisTaskFlowStrange::ReadFromESD(AliESDEvent *tESD) {
1111 AliStack *stack=NULL;
1113 AliMCEvent *mcevent=NULL;
1114 mcevent = MCEvent();
1115 if(mcevent) stack = mcevent->Stack();
1118 Int_t num = tESD->GetNumberOfTracks();
1119 AliESDtrack *myTrack;
1120 Int_t plist[3000], nlist[3000], np=0, nn=0;
1121 Double_t pd0[3000], nd0[3000];
1122 for (Int_t i=0; i!=num; ++i) {
1123 myTrack = (AliESDtrack*) tESD->GetTrack(i);
1124 if(!myTrack) continue;
1126 FillTrackSpy("ESD_TrkAll");
1127 if(!AcceptDaughter()) continue;
1128 FillTrackSpy("ESD_TrkSel");
1129 ((TH2D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("PtIPXY" ))->Fill( myTrack->Pt(), fDaughterImpactParameterXY );
1130 if( myTrack->Charge()>0 ) {
1131 pd0[np] = fDaughterImpactParameterXY;
1134 nd0[nn] = fDaughterImpactParameterXY;
1138 ((TH1D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("NPAIR" ))->Fill( np,nn );
1139 const AliESDVertex *vtx = tESD->GetPrimaryVertex();
1140 AliESDtrack *pT, *nT;
1141 for(int p=0; p!=np; ++p) {
1142 pT = (AliESDtrack*) tESD->GetTrack( plist[p] );
1143 for(int n=0; n!=nn; ++n) {
1144 nT = (AliESDtrack*) tESD->GetTrack( nlist[n] );
1145 fDecayProductIPXY = pd0[p]*nd0[n];
1146 AliExternalTrackParam pETP(*pT), nETP(*nT);
1148 pETP.GetDCA(&nETP,tESD->GetMagneticField(),xa,xb);
1149 fDecayDCAdaughters = pETP.PropagateToDCA(&nETP,tESD->GetMagneticField());
1150 AliESDv0 vertex( nETP,nlist[n], pETP,plist[p] );
1151 fDecayCosinePointingAngleXY = CosThetaPointXY( &vertex, vtx );
1152 fDecayRadXY = DecayLengthXY( &vertex, vtx );
1153 fDecayPt = vertex.Pt();
1154 fDecayPhi = vertex.Phi();
1155 fDecayEta = vertex.Eta();
1156 Double_t pmx, pmy, pmz, nmx, nmy, nmz;
1157 vertex.GetNPxPyPz(nmx,nmy,nmz);
1158 vertex.GetPPxPyPz(pmx,pmy,pmz);
1159 TVector3 mom1(pmx,pmy,pmz), mom2(nmx,nmy,nmz), mom(vertex.Px(),vertex.Py(),vertex.Pz());
1160 Double_t qlpos = mom1.Dot(mom)/mom.Mag();
1161 Double_t qlneg = mom2.Dot(mom)/mom.Mag();
1162 fDecayQt = mom1.Perp(mom);
1163 fDecayAlpha = (qlpos-qlneg)/(qlpos+qlneg);
1164 Double_t mpi = 0.13957018;
1166 Double_t eppi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1167 Double_t enpi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1168 fDecayMass = TMath::Sqrt( mpi*mpi + mpi*mpi + 2*(eppi*enpi - pmx*nmx - pmy*nmy - pmz*nmz ) );
1169 fDecayRapidity = vertex.RapK0Short();
1171 Double_t mpr = 0.938272013;
1174 epr = TMath::Sqrt( mpr*mpr + pmx*pmx + pmy*pmy + pmz*pmz );
1175 epi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1177 epi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1178 epr = TMath::Sqrt( mpr*mpr + nmx*nmx + nmy*nmy + nmz*nmz );
1180 fDecayMass = TMath::Sqrt( mpi*mpi + mpr*mpr + 2*(epi*epr - pmx*nmx - pmy*nmy - pmz*nmz ) );
1181 fDecayRapidity = vertex.RapLambda();
1183 Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + vertex.Px()*vertex.Px() + vertex.Py()*vertex.Py() + vertex.Pz()*vertex.Pz() );
1184 Double_t gamma = energy/fDecayMass;
1185 fDecayDecayLength = DecayLength( &vertex, vtx )/gamma;
1186 fDecayDecayLengthLab = DecayLength( &vertex, vtx );
1187 Double_t dPHI = fDecayPhi;
1188 Double_t dDPHI = dPHI - fPsi2;
1189 if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1190 if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1192 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1193 else FillCandidateSpy("V0SAllIP");
1195 FillCandidateSpy("V0SAll");
1196 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1197 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1198 if(!AcceptCandidate()) continue;
1199 if(fDecayMass<fMinMass) continue;
1200 if(fDecayMass>fMaxMass) continue;
1203 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1204 else FillCandidateSpy("V0SSelIP");
1206 FillCandidateSpy("V0SSel");
1207 ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1208 ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1210 fDecayIDneg = nT->GetID();
1211 fDecayIDpos = pT->GetID();
1212 if(fUseFP) MakeTrack();
1213 LoadTrack(pT); FillTrackSpy("SelDau");
1214 LoadTrack(nT); FillTrackSpy("SelDau");
1216 //===== BEGIN OF MCMATCH
1218 bool matched = false;
1219 Int_t labelpos = pT->GetLabel();
1220 Int_t labelneg = nT->GetLabel();
1222 if( labelpos>0 && labelneg>0 ) {
1223 TParticle *mcpos = stack->Particle( labelpos );
1224 TParticle *mcneg = stack->Particle( labelneg );
1225 Int_t pdgRecPos = mcpos->GetPdgCode();
1226 Int_t pdgRecNeg = mcneg->GetPdgCode();
1227 if( pdgRecPos==211&&pdgRecNeg==-211 ) if(mcpos->GetMother(0)>0) {
1228 if( mcpos->GetMother(0)==mcneg->GetMother(0) ) {
1229 TParticle *mcmot = stack->Particle( mcpos->GetMother(0) );
1230 rOri = TMath::Sqrt( mcmot->Vx()*mcmot->Vx() + mcmot->Vy()*mcmot->Vy() );
1231 if( TMath::Abs(mcmot->GetPdgCode())==310) {
1232 if(mcmot->GetNDaughters()==2) matched=true;
1238 FillCandidateSpy("Mth");
1239 ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1240 ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1241 ((TH1D*)((TList*)fList->FindObject("Mth"))->FindObject("MCOrigin"))->Fill( rOri );
1242 LoadTrack(pT); FillTrackSpy("MthDau");
1243 LoadTrack(nT); FillTrackSpy("MthDau");
1246 //===== END OF MCMATCH
1250 //=======================================================================
1251 void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {
1252 if(!mcArray) return;
1253 AliAODMCParticle *myMCTrack;//, *iMCDau, *jMCDau;
1254 for(int i=0; i!=mcArray->GetEntriesFast(); ++i) {
1255 myMCTrack = dynamic_cast<AliAODMCParticle*>(mcArray->At( i ));
1256 if(!myMCTrack) continue;
1259 if(fSpecie>0) tPDG = 3122;
1260 if( TMath::Abs(myMCTrack->PdgCode())==tPDG )
1261 if( myMCTrack->GetNDaughters() == 2 ) {
1262 Int_t iDau = myMCTrack->GetDaughter(0);
1263 Int_t jDau = myMCTrack->GetDaughter(1);
1264 AliAODMCParticle *posDau=NULL;
1265 AliAODMCParticle *negDau=NULL;
1266 if(iDau>0&&jDau>0) {
1267 iMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( iDau ));
1268 jMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( jDau ));
1270 if(iMCDau->Charge()>0) posDau=iMCDau;
1274 if(jMCDau->Charge()>0) posDau=jMCDau;
1277 } //got two daughters
1278 if(posDau&&negDau) {
1279 Double_t dx = myMCTrack->Xv() - posDau->Xv();
1280 Double_t dy = myMCTrack->Yv() - posDau->Yv();
1281 Double_t dz = myMCTrack->Zv() - posDau->Zv();
1282 fDecayRadXY = TMath::Sqrt( dx*dx + dy*dy );
1283 TVector3 momPos(posDau->Px(),posDau->Py(),posDau->Pz());
1284 TVector3 momNeg(negDau->Px(),negDau->Py(),negDau->Pz());
1285 TVector3 momTot(myMCTrack->Px(),myMCTrack->Py(),myMCTrack->Pz());
1286 Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();
1287 Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();
1288 fDecayQt = momPos.Perp(momTot);
1289 fDecayAlpha = 1.-2./(1.+qlpos/qlneg);
1290 fDecayMass = myMCTrack->GetCalcMass();
1291 Double_t energy = myMCTrack->E();
1292 Double_t gamma = energy/fDecayMass;
1293 fDecayDecayLength = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/gamma;
1294 fDecayPt = myMCTrack->Pt();
1295 fDecayPhi = myMCTrack->Phi();
1296 fDecayEta = myMCTrack->Eta();
1297 fDecayRapidity = myMCTrack->Y();
1298 fDecayDCAdaughters = 0;
1299 fDecayCosinePointingAngleXY = 1;
1300 fDecayProductIPXY = -1;
1301 if(AcceptCandidate()) FillCandidateSpy("GenTru");
1303 } // k0/lda with two daughters
1305 //==== BEGIN TRACK CUTS
1306 if(myMCTrack->Eta()<-0.8) continue;
1307 if(myMCTrack->Eta()>+0.8) continue;
1308 if(myMCTrack->Y()<-0.5) continue;
1309 if(myMCTrack->Y()>+0.5) continue;
1310 //==== END TRACK CUTS
1311 switch( TMath::Abs(myMCTrack->PdgCode()) ) {
1313 FillMCParticleSpy( "MCTPion", myMCTrack );
1314 if( myMCTrack->IsPrimary() )
1315 FillMCParticleSpy( "MCTPionGenAcc", myMCTrack );
1318 FillMCParticleSpy( "MCTKaon", myMCTrack );
1319 if( myMCTrack->IsPrimary() )
1320 FillMCParticleSpy( "MCTKaonGenAcc", myMCTrack );
1323 FillMCParticleSpy( "MCTK0s", myMCTrack );
1324 if( myMCTrack->IsPrimary() )
1325 FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );
1327 case (2212): //proton
1328 FillMCParticleSpy( "MCTProton", myMCTrack );
1329 if( myMCTrack->IsPrimary() )
1330 FillMCParticleSpy( "MCTProtonGenAcc", myMCTrack );
1333 FillMCParticleSpy( "MCTLda", myMCTrack );
1334 if( myMCTrack->IsPrimary() )
1335 FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );
1338 if( myMCTrack->IsPrimary() )
1339 FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );
1342 if( myMCTrack->IsPrimary() )
1343 FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );
1345 case (3334): //omega
1346 if( myMCTrack->IsPrimary() )
1347 FillMCParticleSpy( "MCTOmegaGenAcc", myMCTrack );
1352 //=======================================================================
1353 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliESDv0 *me, const AliVVertex *vtx) {
1354 TVector3 mom( me->Px(), me->Py(), 0 );
1355 TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1356 Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1359 //=======================================================================
1360 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliAODv0 *me, const AliVVertex *vtx) {
1361 TVector3 mom( me->Px(), me->Py(), 0 );
1362 TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1363 Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1366 //=======================================================================
1367 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliESDv0 *me, const AliVVertex *vtx) {
1368 Double_t dx = me->Xv()-vtx->GetX();
1369 Double_t dy = me->Yv()-vtx->GetY();
1370 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1373 //=======================================================================
1374 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliAODv0 *me, const AliVVertex *vtx) {
1375 Double_t dx = me->Xv()-vtx->GetX();
1376 Double_t dy = me->Yv()-vtx->GetY();
1377 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1380 //=======================================================================
1381 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliESDv0 *me, const AliVVertex *vtx) {
1382 Double_t dx = me->Xv()-vtx->GetX();
1383 Double_t dy = me->Yv()-vtx->GetY();
1384 Double_t dz = me->Zv()-vtx->GetZ();
1385 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1388 //=======================================================================
1389 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliAODv0 *me, const AliVVertex *vtx) {
1390 Double_t dx = me->Xv()-vtx->GetX();
1391 Double_t dy = me->Yv()-vtx->GetY();
1392 Double_t dz = me->Zv()-vtx->GetZ();
1393 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1396 //=======================================================================
1397 void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
1398 TClonesArray* mcArray=NULL;
1400 mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1404 Int_t nV0s = tAOD->GetNumberOfV0s();
1406 Int_t v0all=0, v0imw=0;
1407 for (Int_t i=0; i!=nV0s; ++i) {
1408 myV0 = (AliAODv0*) tAOD->GetV0(i);
1410 if(!fOnline) if(myV0->GetOnFlyStatus() ) continue;
1411 if(fOnline) if(!myV0->GetOnFlyStatus() ) continue;
1413 fDecayPt = myV0->Pt();
1414 fDecayPhi = myV0->Phi();
1415 fDecayEta = myV0->Eta();
1417 AliAODTrack *iT, *jT;
1418 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1419 Double_t pos[3],cov[6];
1421 vtx->GetCovarianceMatrix(cov);
1422 const AliESDVertex vESD(pos,cov,100.,100);
1425 iT=(AliAODTrack*) myV0->GetDaughter(0);
1426 if(iT->Charge()>0) {
1433 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1434 AliESDtrack ieT( iT );
1435 ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1436 ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1437 ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1438 ieT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1439 LoadTrack(&ieT,iT->Chi2perNDF());
1441 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1442 fDaughterImpactParameterXY = ip[0];
1443 fDaughterImpactParameterZ = ip[1];
1444 fDecayIPpos = fDaughterImpactParameterXY; //ieT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1445 FillTrackSpy("AllDau");
1446 if(!AcceptDaughter(fDecayPt<2.0?kTRUE:kFALSE)) continue;
1448 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1449 AliESDtrack jeT( jT );
1450 jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1451 jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1452 jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1453 jeT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1454 LoadTrack(&jeT,jT->Chi2perNDF());
1455 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1456 fDaughterImpactParameterXY = ip[0];
1457 fDaughterImpactParameterZ = ip[1];
1458 fDecayIPneg = fDaughterImpactParameterXY; //jeT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1459 FillTrackSpy("AllDau");
1460 if(!AcceptDaughter(fDecayPt<2.0?kTRUE:kFALSE)) continue;
1462 if( fExcludeTPCEdges ) {
1463 if( IsAtTPCEdge(iT->Phi(),iT->Pt(),+1,tAOD->GetMagneticField()) ) continue;
1464 if( IsAtTPCEdge(jT->Phi(),jT->Pt(),-1,tAOD->GetMagneticField()) ) continue;
1466 ieT.GetDCA(&jeT,tAOD->GetMagneticField(),fDecayXpos,fDecayXneg);
1468 // cutting out population close to TPC edges :: strange excess saw in 2010
1469 if( fExcludeTPCEdges ) {
1470 Double_t phimod = myV0->Phi();
1471 int sectors[6] = {5,6,9,10,11,12};
1472 for(int ii=0; ii!=6; ++ii)
1473 if( (phimod<(sectors[ii]+1)*TMath::Pi()/9.0) && (phimod>sectors[ii]*TMath::Pi()/9.0) )
1478 fDecayRapidity = myV0->RapK0Short();
1480 fDecayRapidity = myV0->RapLambda();
1481 fDecayDCAdaughters = myV0->DcaV0Daughters();
1482 fDecayCosinePointingAngleXY = CosThetaPointXY( myV0, vtx );
1483 fDecayRadXY = DecayLengthXY( myV0, vtx );
1484 fDecayProductIPXY = fDecayIPpos*fDecayIPneg;
1485 fDecayQt = myV0->PtArmV0();
1486 fDecayAlpha = myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1487 if(myV0->ChargeProng(iPos)<0) fDecayAlpha = -fDecayAlpha; // protects for a change in convention
1488 fDecayPt = myV0->Pt();
1489 fDecayEta = myV0->Eta();
1491 fDecayMass = myV0->MassK0Short();
1493 if(fDecayAlpha>0) fDecayMass = myV0->MassLambda();
1494 else fDecayMass = myV0->MassAntiLambda();
1497 if(fDecayMass<fMinMass) continue;
1498 if(fDecayMass>fMaxMass) continue;
1500 Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + myV0->Px()*myV0->Px() + myV0->Py()*myV0->Py() + myV0->Pz()*myV0->Pz() );
1501 Double_t gamma = energy/fDecayMass;
1502 fDecayDecayLength = DecayLength( myV0, vtx )/gamma;
1503 fDecayDecayLengthLab = DecayLength( myV0, vtx );
1504 Double_t dPHI = fDecayPhi;
1505 Double_t dDPHI = dPHI - fPsi2;
1506 if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1507 if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1509 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1510 else FillCandidateSpy("V0SAllIP");
1512 FillCandidateSpy("V0SAll");
1514 FillDecayVn("V0SAllVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1516 if(!AcceptCandidate()) continue;
1518 if(fDecayPt<fDecayStopPIDAtPt) {
1519 if( fSpecie==0 ) {//PID for kzero::pion+pion
1520 if( !PassesPIDCuts(&ieT,AliPID::kPion) ) continue; //positive track
1521 if( !PassesPIDCuts(&jeT,AliPID::kPion) ) continue; //negative track
1522 } else { //PID for lambda::proton+pion
1524 if( !PassesPIDCuts(&ieT,AliPID::kProton) ) continue; //positive track
1525 if( !PassesPIDCuts(&jeT,AliPID::kPion) ) continue; //negative track
1527 if( !PassesPIDCuts(&jeT,AliPID::kProton) ) continue; //negative track
1528 if( !PassesPIDCuts(&ieT,AliPID::kPion) ) continue; //positive track
1533 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1534 else FillCandidateSpy("V0SSelIP");
1536 FillCandidateSpy("V0SSel");
1538 FillDecayVn("V0SSelVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1539 // ============================
1540 // Posting for FlowAnalysis
1542 fDecayIDneg = iT->GetID();
1543 fDecayIDpos = jT->GetID();
1544 if(fUseFP) MakeTrack();
1546 // ============================
1547 LoadTrack(&ieT,iT->Chi2perNDF());
1548 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1549 fDaughterImpactParameterXY = ip[0];
1550 fDaughterImpactParameterZ = ip[1];
1551 FillTrackSpy("SelDau");
1552 LoadTrack(&jeT,jT->Chi2perNDF());
1553 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1554 fDaughterImpactParameterXY = ip[0];
1555 fDaughterImpactParameterZ = ip[1];
1556 FillTrackSpy("SelDau");
1557 //===== BEGIN OF MCMATCH
1558 if(fReadMC) ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event
1560 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
1561 bool matched = false;
1562 bool feeddown = false;
1563 Int_t labelpos = iT->GetLabel();
1564 Int_t labelneg = jT->GetLabel();
1565 AliAODMCParticle *mcpos = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelpos) );
1566 AliAODMCParticle *mcneg = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelneg) );
1567 if( mcpos && mcneg ) {
1568 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Daughters in stack
1569 Int_t pdgRecPos = mcpos->GetPdgCode();
1570 Int_t pdgRecNeg = mcneg->GetPdgCode();
1571 int pospdg=211, negpdg=211;
1572 int mompdg=310, fdwpdg=333;
1577 pospdg=2212; negpdg=211;
1579 negpdg=2212; pospdg=211;
1582 if( TMath::Abs(pdgRecPos)==pospdg&&TMath::Abs(pdgRecNeg)==negpdg )
1583 if(mcpos->GetMother()>-1)
1584 if( mcpos->GetMother()==mcneg->GetMother() ) {
1585 AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );
1586 fDecayMatchOrigin = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );
1587 fDecayMatchPt = mcmot->Pt();
1588 fDecayMatchEta = mcmot->Eta();
1589 fDecayMatchPhi = mcmot->Phi();
1590 if( TMath::Abs(mcmot->GetPdgCode())==mompdg) {
1591 if(mcmot->GetNDaughters()==2) {
1592 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 4 ); // Correspond to decay
1594 Double_t dx = mcmot->Xv() - mcpos->Xv();
1595 Double_t dy = mcmot->Yv() - mcpos->Yv();
1596 fDecayMatchRadXY = TMath::Sqrt( dx*dx + dy*dy );
1598 if(mcmot->GetMother()>-1) {
1599 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 5 ); // Decay has mother
1600 AliAODMCParticle *mcfdw = (AliAODMCParticle*) mcArray->At( mcmot->GetMother() );
1601 if( TMath::Abs(mcfdw->GetPdgCode())==fdwpdg)
1603 } // k0/lda have mother
1604 } // mother matches k0/lda
1605 } // both have same mother
1608 FillCandidateSpy("Mth",true);
1610 FillDecayVn("V0SMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1611 if(fPostMatched>0) {
1612 fDecayIDneg = iT->GetID();
1613 fDecayIDpos = jT->GetID();
1614 if(fUseFP) MakeTrack();
1616 if(labelpos<0&&labelneg<0) {
1617 FillCandidateSpy("MthNegNeg",true);
1619 FillDecayVn("V0SMthNegNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1620 } else if(labelpos>0&&labelneg>0) {
1622 FillDecayVn("V0SMthPosPosVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1623 } else if(labelpos*labelneg<0) {
1624 FillCandidateSpy("MthPosNeg",true);
1626 FillDecayVn("V0SMthPosNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1628 AliAODVertex *secvtx = myV0->GetSecondaryVtx();
1629 Double_t possec[3],covsec[6];
1630 secvtx->GetXYZ(possec);
1631 secvtx->GetCovarianceMatrix(covsec);
1632 const AliESDVertex vSecVtx(possec,covsec,100.,100);
1633 AliESDtrack trackAtSecI( iT );
1634 trackAtSecI.SetTPCClusterMap( iT->GetTPCClusterMap() );
1635 trackAtSecI.SetTPCSharedMap( iT->GetTPCSharedMap() );
1636 trackAtSecI.SetTPCPointsF( iT->GetTPCNclsF() );
1637 trackAtSecI.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1638 fDaughterAtSecPhi = trackAtSecI.Phi();
1639 fDaughterAtSecEta = trackAtSecI.Eta();
1640 fDaughterAtSecPt = trackAtSecI.Pt();
1641 LoadTrack(&ieT,iT->Chi2perNDF());
1642 fDaughterMatchPhi=mcpos->Phi();
1643 fDaughterMatchEta=mcpos->Eta();
1644 fDaughterMatchPt=mcpos->Pt();
1645 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1646 fDaughterImpactParameterXY = ip[0];
1647 fDaughterImpactParameterZ = ip[1];
1648 FillTrackSpy("MthDau",true);
1649 if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1650 else FillTrackSpy("MthPosDau",true);
1651 AliESDtrack trackAtSecJ( jT );
1652 trackAtSecJ.SetTPCClusterMap( jT->GetTPCClusterMap() );
1653 trackAtSecJ.SetTPCSharedMap( jT->GetTPCSharedMap() );
1654 trackAtSecJ.SetTPCPointsF( jT->GetTPCNclsF() );
1655 trackAtSecJ.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1656 fDaughterAtSecPhi = trackAtSecJ.Phi();
1657 fDaughterAtSecEta = trackAtSecJ.Eta();
1658 fDaughterAtSecPt = trackAtSecJ.Pt();
1659 LoadTrack(&jeT,jT->Chi2perNDF());
1660 fDaughterMatchPhi=mcneg->Phi();
1661 fDaughterMatchEta=mcneg->Eta();
1662 fDaughterMatchPt=mcneg->Pt();
1663 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1664 fDaughterImpactParameterXY = ip[0];
1665 fDaughterImpactParameterZ = ip[1];
1666 FillTrackSpy("MthDau",true);
1667 if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1668 else FillTrackSpy("MthPosDau",true);
1670 FillCandidateSpy("UnMth",false);
1672 FillDecayVn("V0SUnMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1673 if(fPostMatched<0) {
1674 fDecayIDneg = iT->GetID();
1675 fDecayIDpos = jT->GetID();
1676 if(fUseFP) MakeTrack();
1678 LoadTrack(&ieT,iT->Chi2perNDF());
1679 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1680 fDaughterImpactParameterXY = ip[0];
1681 fDaughterImpactParameterZ = ip[1];
1682 FillTrackSpy("UnMthDau",false);
1683 LoadTrack(&jeT,jT->Chi2perNDF());
1684 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1685 fDaughterImpactParameterXY = ip[0];
1686 fDaughterImpactParameterZ = ip[1];
1687 FillTrackSpy("UnMthDau",false);
1690 FillCandidateSpy("MthFeedDown",true);
1693 //===== END OF MCMATCH
1695 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("V0SADC"))->Fill( v0all,v0imw );
1697 QCStoreDecayVn("V0SAllVn");
1698 QCStoreDecayVn("V0SSelVn");
1700 QCStoreDecayVn("V0SMthVn");
1701 QCStoreDecayVn("V0SMthNegNegVn");
1702 QCStoreDecayVn("V0SMthPosPosVn");
1703 QCStoreDecayVn("V0SMthPosNegVn");
1708 //=======================================================================
1709 Bool_t AliAnalysisTaskFlowStrange::PassesPIDCuts(AliESDtrack *myTrack, AliPID::EParticleType pid) {
1712 fDaughterNSigmaPID = fPIDResponse->NumberOfSigmasTPC(myTrack,pid);
1713 if( TMath::Abs(fDaughterNSigmaPID) > fDaughterMaxNSigmaPID )
1718 //=======================================================================
1719 void AliAnalysisTaskFlowStrange::ChargeParticles(AliAODEvent *tAOD) {
1720 //benchmark purposes
1722 TClonesArray* mcArray=NULL;
1724 mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1727 for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
1728 AliAODTrack *t = dynamic_cast<AliAODTrack*>(tAOD->GetTrack( i ));
1730 if( !t->TestFilterBit(1) ) continue;
1731 fDecayMass=0.0; // using mass as nsigmas control plot
1732 if(fPIDResponse) { // PID
1733 switch(fSpecie) { // TPC PID only
1735 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kPion);
1738 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kKaon);
1741 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kProton);
1745 Bool_t pass = kTRUE;
1746 if( TMath::Abs(fDecayMass) > 3.0 ) pass=kFALSE;
1747 if( t->Eta()<-0.5 || t->Eta()>+0.5 ) pass=kFALSE;
1748 if( t->Pt()<0.2 || t->Pt()>20.0 ) pass=kFALSE;
1749 AliESDtrack et( t );
1750 et.SetTPCClusterMap( t->GetTPCClusterMap() );
1751 et.SetTPCSharedMap( t->GetTPCSharedMap() );
1752 et.SetTPCPointsF( t->GetTPCNclsF() );
1754 LoadTrack(&et,t->Chi2perNDF());
1755 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1758 et.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1759 fDaughterImpactParameterXY = ip[0];
1760 fDaughterImpactParameterZ = ip[1];
1762 FillTrackSpy("TrkAll");
1764 FillTrackVn("TrkAllVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1766 FillTrackSpy("TrkSel");
1768 FillTrackVn("TrkSelVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1770 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event
1772 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
1773 bool matched = false;
1774 Int_t label = t->GetLabel();
1775 AliAODMCParticle *mcpar = (AliAODMCParticle*) mcArray->At( TMath::Abs(label) );
1777 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Particle in stack
1778 Int_t pdgmcpar = TMath::Abs(mcpar->GetPdgCode());
1781 if(pdgmcpar==211) matched = true;
1784 if(pdgmcpar==211) matched = true;
1787 if(pdgmcpar==2212) matched = true;
1790 if(!mcpar->IsPrimary()) matched = false;
1793 FillTrackSpy("Mth");
1795 FillTrackVn("MthVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1797 FillTrackSpy("MthNeg");
1799 FillTrackVn("MthNegVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1801 FillTrackSpy("MthPos");
1803 FillTrackVn("MthPosVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1812 fDecayID=t->GetID();
1817 QCStoreTrackVn("TrkAllVn");
1818 QCStoreTrackVn("TrkSelVn");
1820 QCStoreTrackVn("MthVn");
1821 QCStoreTrackVn("MthNegVn");
1822 QCStoreTrackVn("MthPosVn");
1827 //=======================================================================
1828 void AliAnalysisTaskFlowStrange::ComputeChi2VZERO() {
1829 Double_t MeanQaQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->GetBinContent( 1 );
1830 Double_t MeanQaQa = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->GetBinContent( 1 );
1831 Double_t MeanQcQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->GetBinContent( 1 );
1832 Double_t MeanQaQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->GetBinContent( 1 );
1833 Double_t MeanQcQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->GetBinContent( 1 );
1834 if(!TMath::AreEqualAbs(MeanQaQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQcQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQaQc,0,1e-10)) {
1835 Double_t OneOverChiSquaredVZEA = MeanQaQa*MeanQcQt/MeanQaQc/MeanQaQt-1;
1836 Double_t OneOverChiSquaredVZEC = MeanQcQc*MeanQaQt/MeanQaQc/MeanQcQt-1;
1837 if(!TMath::AreEqualAbs(OneOverChiSquaredVZEA,0,1e-10)&&!TMath::AreEqualAbs(OneOverChiSquaredVZEC,0,1e-10)) {
1838 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->SetBinContent( 1, 1/OneOverChiSquaredVZEA );
1839 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->SetBinContent( 1, 1/OneOverChiSquaredVZEC );
1843 //=======================================================================
1844 void AliAnalysisTaskFlowStrange::Terminate(Option_t *) {
1846 if(fSkipTerminate) return;
1848 if(fSkipSelection) return;
1851 ComputeDecayVn("V0SAllVn");
1852 ComputeDecayVn("V0SSelVn");
1854 ComputeDecayVn("V0SMthVn");
1855 ComputeDecayVn("V0SMthPosPosVn");
1856 ComputeDecayVn("V0SMthNegNegVn");
1857 ComputeDecayVn("V0SMthPosNegVn");
1858 ComputeDecayVn("V0SUnMthVn");
1861 ComputeTrackVn("TrkAllVn");
1862 ComputeTrackVn("TrkSelVn");
1864 ComputeTrackVn("MthVn");
1865 ComputeTrackVn("MthPosVn");
1866 ComputeTrackVn("MthNegVn");
1870 //=======================================================================
1871 void AliAnalysisTaskFlowStrange::MakeTrack() {
1872 // create track for flow tasks
1873 if(fCandidates->GetLast()+5>fCandidates->GetSize()) {
1874 fCandidates->Expand( fCandidates->GetSize()+20 );
1876 Bool_t overwrite = kTRUE;
1877 AliFlowCandidateTrack *oTrack = (static_cast<AliFlowCandidateTrack*> (fCandidates->At( fCandidates->GetLast()+1 )));
1878 if( !oTrack ) { // creates new
1879 oTrack = new AliFlowCandidateTrack();
1881 } else { // overwrites
1884 oTrack->SetMass(fDecayMass);
1885 oTrack->SetPt(fDecayPt);
1886 oTrack->SetPhi(fDecayPhi);
1887 oTrack->SetEta(fDecayEta);
1889 oTrack->AddDaughter(fDecayIDpos);
1890 oTrack->AddDaughter(fDecayIDneg);
1892 oTrack->SetID( fDecayID );
1894 oTrack->SetForPOISelection(kTRUE);
1895 oTrack->SetForRPSelection(kFALSE);
1897 fCandidates->SetLast( fCandidates->GetLast()+1 );
1899 fCandidates->AddLast(oTrack);
1903 //=======================================================================
1904 void AliAnalysisTaskFlowStrange::AddCandidates() {
1905 // adds candidates to flow events (untaging if necessary)
1906 if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1907 if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1908 if(fDebug) printf("I received %d candidates\n",fCandidates->GetEntriesFast());
1911 for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
1912 AliFlowCandidateTrack *cand = static_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
1914 cand->SetForPOISelection(kTRUE);
1915 cand->SetForRPSelection(kFALSE);
1917 //if(fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
1918 if(fSpecie<10) { // DECAYS
1920 if(fDaughterUnTag) {
1921 for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
1922 if(fDebug) printf(" >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
1923 for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1924 AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1926 if(!iRP->InRPSelection()) continue;
1927 if(cand->GetIDDaughter(iDau) == iRP->GetID()) {
1928 if(fDebug) printf(" was in RP set");
1930 iRP->SetForRPSelection(kFALSE);
1931 fTPCevent->SetNumberOfRPs( fTPCevent->GetNumberOfRPs() -1 );
1934 if(fDebug) printf("\n");
1938 fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1940 // adding only new tracks and tagging accordingly ===>
1941 Bool_t found=kFALSE;
1942 for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1943 AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1945 if(!iRP->InRPSelection()) continue;
1946 if(cand->GetID() == iRP->GetID()) {
1947 if(fDebug) printf(" >charged track (%d) was also found in RP set (adding poi tag)\n",cand->GetID());
1948 iRP->SetMass( cand->Mass() );
1949 iRP->SetForPOISelection(kTRUE);
1953 if(!found) // not found adding track
1954 fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1956 fVZEevent->InsertTrack( ((AliFlowTrack*) cand) );
1958 fTPCevent->SetNumberOfPOIs( poi );
1959 fVZEevent->SetNumberOfPOIs( poi );
1960 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("POI"))->Fill( poi );
1961 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("UNTAG"))->Fill( untagged );
1962 if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1963 if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1965 //=======================================================================
1966 void AliAnalysisTaskFlowStrange::PushBackFlowTrack(AliFlowEvent *flowevent, Double_t pt, Double_t phi, Double_t eta, Double_t w, Int_t id) {
1972 rfp.SetForRPSelection(kTRUE);
1973 rfp.SetForPOISelection(kFALSE);
1976 flowevent->InsertTrack( &rfp );
1978 //=======================================================================
1979 Bool_t AliAnalysisTaskFlowStrange::IsAtTPCEdge(Double_t phi,Double_t pt,Int_t charge,Double_t b) {
1980 // Origin: Alex Dobrin
1981 // Implemented by Carlos Perez
1982 TF1 cutLo("cutLo", "-0.01/x+pi/18.0-0.015", 0, 100);
1983 TF1 cutHi("cutHi", "0.55/x/x+pi/18.0+0.03", 0, 100);
1984 Double_t phimod = phi;
1985 if(b<0) phimod = TMath::TwoPi()-phimod; //for negatve polarity field
1986 if(charge<0) phimod = TMath::TwoPi()-phimod; //for negatve charge
1987 phimod += TMath::Pi()/18.0;
1988 phimod = fmod(phimod, TMath::Pi()/9.0);
1989 if( phimod<cutHi.Eval(pt) && phimod>cutLo.Eval(pt) )
1994 //=======================================================================
1995 void AliAnalysisTaskFlowStrange::MakeQVectors() {
1996 //computes event plane and updates fPsi2
1997 //if there is a problem fPsi->-1
2001 MakeQVZE(InputEvent());
2002 MakeQTPC(InputEvent());
2003 if(fUseFP&&fReadMC) {
2004 fVZEevent->SetMCReactionPlaneAngle( fMCEP );
2005 fTPCevent->SetMCReactionPlaneAngle( fMCEP );
2008 printf("**::MakeQVectors()");
2009 printf(" fQVZEACos %.16f | fQVZEASin %.16f || fQVZEA %.3f | fQVZEC %.3f \n",fQVZEACos, fQVZEASin, fQVZEA, fQVZEC);
2010 printf(" nQTPA_nTracks %d | fQTPC_nTracks %d || fQTPCA %.3f | fQTPCC %.3f \n",fQTPCA_nTracks, fQTPCC_nTracks, fQTPCA, fQTPCC);
2011 printf(" fQTPCACos %.16f | fQTPCASin %.16f || fQTPC2hCos %.16f | fQTPC2hSin %.16f \n",fQTPCACos, fQTPCASin, fQTPC2hCos, fQTPC2hSin);
2015 //=======================================================================
2016 void AliAnalysisTaskFlowStrange::FillMakeQSpy() {
2019 Double_t qvzecos,qvzesin,psivzea,psivzec,psivze,qvze, qvzea, qvzec;
2020 psivzea = ( TMath::Pi()+TMath::ATan2(-fQVZEASin,-fQVZEACos) )/fHarmonic;
2021 psivzec = ( TMath::Pi()+TMath::ATan2(-fQVZECSin,-fQVZECCos) )/fHarmonic;
2022 qvzecos = fQVZEACos + fQVZECCos;
2023 qvzesin = fQVZEASin + fQVZECSin;
2026 qvze = fQVZEA + fQVZEC;
2027 psivze = ( TMath::Pi()+TMath::ATan2(-qvzesin,-qvzecos) )/fHarmonic;
2029 Double_t qtpccos,qtpcsin,psitpca,psitpcc,psitpc,qtpc;
2030 psitpca = ( TMath::Pi()+TMath::ATan2(-fQTPCASin,-fQTPCACos) )/fHarmonic;
2031 psitpcc = ( TMath::Pi()+TMath::ATan2(-fQTPCCSin,-fQTPCCCos) )/fHarmonic;
2032 qtpccos = fQTPCACos + fQTPCCCos;
2033 qtpcsin = fQTPCASin + fQTPCCSin;
2034 qtpc = fQTPCA + fQTPCC;
2035 psitpc = ( TMath::Pi()+TMath::ATan2(-qtpcsin,-qtpccos) )/fHarmonic;
2036 //=>does the event clear?
2039 if(fVZEWarning) return;
2043 if(fQTPCA<2||fQTPCC<2) return;
2047 //computing physical Qm vectors
2048 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
2049 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
2050 Double_t vzec_qmnor = TMath::Sqrt( vzec_qmcos*vzec_qmcos + vzec_qmsin*vzec_qmsin );
2051 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
2052 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
2053 Double_t vzea_qmnor = TMath::Sqrt( vzea_qmcos*vzea_qmcos + vzea_qmsin*vzea_qmsin );
2054 Double_t vze_qmcos = qvzecos/qvze;
2055 Double_t vze_qmsin = qvzesin/qvze;
2056 Double_t vze_qmnor = TMath::Sqrt( vze_qmcos*vze_qmcos + vze_qmsin*vze_qmsin );
2057 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
2058 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
2059 Double_t tpcc_qmnor = TMath::Sqrt( tpcc_qmcos*tpcc_qmcos + tpcc_qmsin*tpcc_qmsin );
2060 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
2061 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
2062 Double_t tpca_qmnor = TMath::Sqrt( tpca_qmcos*tpca_qmcos + tpca_qmsin*tpca_qmsin );
2063 Double_t tpc_qmcos = qtpccos/qtpc;
2064 Double_t tpc_qmsin = qtpcsin/qtpc;
2065 Double_t tpc_qmnor = TMath::Sqrt( tpc_qmcos*tpc_qmcos + tpc_qmsin*tpc_qmsin );
2066 //=>great! recording
2067 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSI"))->Fill( psivze );
2068 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIA"))->Fill( psivzea );
2069 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIC"))->Fill( psivzec );
2070 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPVZE"))->Fill( qvze );
2071 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEA"))->Fill( vzea_qmnor*TMath::Sqrt(qvzea) );
2072 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEC"))->Fill( vzec_qmnor*TMath::Sqrt(qvzec) );
2074 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSI"))->Fill( psitpc );
2075 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIA"))->Fill( psitpca );
2076 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIC"))->Fill( psitpcc );
2077 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPTPC"))->Fill( qtpc );
2078 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPC"))->Fill( tpc_qmnor*TMath::Sqrt(qtpc) );
2080 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCAVZEC"))->Fill( psitpca, psivzec );
2081 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCCVZEA"))->Fill( psitpcc, psivzea );
2082 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCVZE"))->Fill( psitpc, psivze );
2085 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPC"))->Fill( psitpc-fMCEP );
2086 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCA"))->Fill( psitpca-fMCEP );
2087 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCC"))->Fill( psitpcc-fMCEP );
2088 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZE"))->Fill( psivze-fMCEP );
2089 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEA"))->Fill( psivzea-fMCEP );
2090 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEC"))->Fill( psivzec-fMCEP );
2093 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 1., tpcc_qmsin, tpcc_qmnor );
2094 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 2., tpcc_qmcos, tpcc_qmnor );
2095 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 3., tpca_qmsin, tpca_qmnor );
2096 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 4., tpca_qmcos, tpca_qmnor );
2097 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 5., tpc_qmsin, tpc_qmnor );
2098 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 6., tpc_qmcos, tpc_qmnor );
2100 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 1., vzec_qmsin, vzec_qmnor );
2101 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 2., vzec_qmcos, vzec_qmnor );
2102 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 3., vzea_qmsin, vzea_qmnor );
2103 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 4., vzea_qmcos, vzea_qmnor );
2104 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 5., vze_qmsin, vze_qmnor );
2105 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 6., vze_qmcos, vze_qmnor );
2107 Double_t vzeqaqc = vzec_qmcos*vzea_qmcos + vzec_qmsin*vzea_qmsin;
2108 Double_t vzeqatpcq = vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin;
2109 Double_t vzeqctpcq = vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin;
2110 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->Fill( 1., vzeqaqc, vze_qmnor );
2111 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->Fill( 1., vzea_qmnor*vzea_qmnor, vze_qmnor );
2112 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->Fill( 1., vzec_qmnor*vzec_qmnor, vze_qmnor );
2113 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->Fill( 1., vzeqatpcq, vze_qmnor );
2114 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->Fill( 1., vzeqctpcq, vze_qmnor );
2117 //=======================================================================
2118 void AliAnalysisTaskFlowStrange::MakeQVZE(AliVEvent *tevent) {
2120 if(fUseFP) fVZEevent->ClearFast(); // flowpackage
2122 fQVZEACos=fQVZEASin=fQVZEA=fQVZECCos=fQVZECSin=fQVZEC=0;
2124 Double_t eta, phi, w;
2126 for(int id=fVZECa*8;id!=8+fVZECb*8;++id) {
2127 eta = -3.45+0.5*(id/8);
2128 phi = TMath::PiOver4()*(0.5+id%8);
2129 w = tevent->GetVZEROEqMultiplicity(id);
2130 if(w<3) fVZEWarning=kTRUE;
2132 fQVZECCos += w*TMath::Cos(fHarmonic*phi);
2133 fQVZECSin += w*TMath::Sin(fHarmonic*phi);
2135 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
2137 if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
2140 for(int id=32+fVZEAa*8;id!=40+fVZEAb*8;++id) {
2141 eta = +4.8-0.6*((id/8)-4);
2142 phi = TMath::PiOver4()*(0.5+id%8);
2143 w = tevent->GetVZEROEqMultiplicity(id);
2144 if(w<3) fVZEWarning=kTRUE;
2146 fQVZEACos += w*TMath::Cos(fHarmonic*phi);
2147 fQVZEASin += w*TMath::Sin(fHarmonic*phi);
2149 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
2151 if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
2153 if(fUseFP) { // flowpackage
2154 fVZEevent->SetNumberOfRPs(rfp);
2155 if(fDebug>0) printf("Inserted tracks in FlowEventVZE %d ==> %.1f\n",fVZEevent->NumberOfTracks(),fQVZEC+fQVZEA);
2158 //=======================================================================
2159 void AliAnalysisTaskFlowStrange::AddTPCRFPSpy(TList *me) {
2161 tH1D = new TH1D("PT", "PT", 50,0,5); me->Add(tH1D);
2162 tH1D = new TH1D("PHI", "PHI", 90,0,TMath::TwoPi()); me->Add(tH1D);
2163 tH1D = new TH1D("ETA", "ETA", 40,-1,+1); me->Add(tH1D);
2164 tH1D = new TH1D("TPCS", "TPC Signal", 100,0,500); me->Add(tH1D);
2165 tH1D = new TH1D("IPXY", "IPXY", 100,-2,+2); me->Add(tH1D);
2166 tH1D = new TH1D("IPZ", "IPZ", 100,-2,+2); me->Add(tH1D);
2168 tH1D = new TH1D("TPCNCLS", "NCLS", 170,-0.5,+169.5); me->Add(tH1D);
2169 tH1D = new TH1D("TPCSHCL", "NSCLS / NCLS", 100,0,1); me->Add(tH1D);
2170 tH1D = new TH1D("TPCFICL", "NCLS1I / NCLS",100,0,1); me->Add(tH1D);
2171 tH1D = new TH1D("TPCXRNF", "XROW / NFCLS", 100,0,1.5); me->Add(tH1D);
2172 tH1D = new TH1D("TPCRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
2174 tH1D = new TH1D("ITSNCLS", "NCLS", 7,-0.5,+6.5); me->Add(tH1D);
2175 tH1D = new TH1D("ITSRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
2178 //=======================================================================
2179 Bool_t AliAnalysisTaskFlowStrange::PassesRFPTPCCuts(AliESDtrack *track, Double_t aodchi2cls, Float_t aodipxy, Float_t aodipz) {
2180 if(track->GetKinkIndex(0)>0) return kFALSE;
2181 if( (track->GetStatus()&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2182 Double_t pt = track->Pt();
2183 Double_t phi = track->Phi();
2184 Double_t eta = track->Eta();
2185 Double_t tpcs = track->GetTPCsignal();
2187 track->GetImpactParameters(ipxy,ipz);
2188 Int_t cls = track->GetTPCclusters(0);
2189 Double_t xrows, findcls, chi2;
2190 findcls = track->GetTPCNclsF();
2191 xrows = track->GetTPCCrossedRows();
2192 chi2 = track->GetTPCchi2();
2193 Double_t rchi2 = chi2/cls;
2199 Double_t xrnfcls = xrows/findcls;
2200 Double_t scls, cls1i, itschi2;
2202 cls1i = track->GetTPCNclsIter1();
2203 scls = track->GetTPCnclsS();
2204 itscls = track->GetITSclusters(0);
2205 itschi2 = track->GetITSchi2();
2206 Double_t shcl = scls/cls;
2207 Double_t ficl = cls1i/cls;
2208 Double_t itsrchi2 = itscls/itschi2;
2209 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PT"))->Fill( pt );
2210 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PHI"))->Fill( phi );
2211 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ETA"))->Fill( eta );
2212 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCS"))->Fill( tpcs );
2213 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPXY"))->Fill( ipxy );
2214 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPZ"))->Fill( ipz );
2215 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCNCLS"))->Fill( cls );
2216 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCSHCL"))->Fill( shcl );
2217 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCFICL"))->Fill( ficl );
2218 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2219 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2220 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSNCLS"))->Fill( itscls );
2221 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2222 if(pt<fRFPminPt) return kFALSE; //0.2
2223 if(pt>fRFPmaxPt) return kFALSE; //5.0
2224 if(eta<fRFPCminEta||(eta>fRFPCmaxEta&&eta<fRFPAminEta)||eta>fRFPAmaxEta) return kFALSE; // -0.8 0.0 0.0 +0.8
2225 if(tpcs<fRFPTPCsignal) return kFALSE; //10.0
2226 if( TMath::Sqrt(ipxy*ipxy/fRFPmaxIPxy/fRFPmaxIPxy+ipz*ipz/fRFPmaxIPz/fRFPmaxIPz)>1 ) return kFALSE; //2.4 3.2
2227 if(cls<fRFPTPCncls) return kFALSE; //70
2228 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PT"))->Fill( pt );
2229 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PHI"))->Fill( phi );
2230 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ETA"))->Fill( eta );
2231 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCS"))->Fill( tpcs );
2232 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPXY"))->Fill( ipxy );
2233 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPZ"))->Fill( ipz );
2234 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCNCLS"))->Fill( cls );
2235 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCSHCL"))->Fill( shcl );
2236 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCFICL"))->Fill( ficl );
2237 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2238 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2239 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSNCLS"))->Fill( itscls );
2240 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2243 //=======================================================================
2244 void AliAnalysisTaskFlowStrange::MakeQTPC(AliVEvent *tevent) {
2245 AliESDEvent *tESD = (AliESDEvent*) (tevent);
2246 AliAODEvent *tAOD = (AliAODEvent*) (tevent);
2255 //=======================================================================
2256 void AliAnalysisTaskFlowStrange::MakeQTPC(AliAODEvent *tAOD) {
2258 if(fUseFP) fTPCevent->ClearFast(); // flowpackage
2259 fQTPCACos=fQTPCASin=fQTPCA=fQTPC2hCos=0;
2260 fQTPCCCos=fQTPCCSin=fQTPCC=fQTPC2hSin=0;
2264 Double_t eta, phi, w;
2266 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
2267 Double_t pos[3],cov[6];
2269 vtx->GetCovarianceMatrix(cov);
2270 const AliESDVertex vESD(pos,cov,100.,100);
2273 Int_t rawN = tAOD->GetNumberOfTracks();
2274 for(Int_t id=0; id!=rawN; ++id) {
2275 track = dynamic_cast<AliAODTrack*>(tAOD->GetTrack(id));
2276 if(!track) AliFatal("Not a standard AOD");
2278 if(!track->TestFilterBit(fRFPFilterBit)) continue;
2279 if( fExcludeTPCEdges )
2280 if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tAOD->GetMagneticField() ) ) continue;
2281 AliESDtrack etrack( track );
2282 etrack.SetTPCClusterMap( track->GetTPCClusterMap() );
2283 etrack.SetTPCSharedMap( track->GetTPCSharedMap() );
2284 etrack.SetTPCPointsF( track->GetTPCNclsF() );
2286 etrack.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
2287 if(!PassesRFPTPCCuts(&etrack,track->Chi2perNDF(),ip[0],ip[1])) continue;
2293 fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2294 fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2296 fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
2298 fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2299 fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2301 fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
2303 fQTPC2hCos += w*TMath::Cos(2.0*fHarmonic*phi);
2304 fQTPC2hSin += w*TMath::Sin(2.0*fHarmonic*phi);
2306 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2307 if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flow package
2310 fTPCevent->SetNumberOfRPs(rfp);
2311 if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
2314 //=======================================================================
2315 void AliAnalysisTaskFlowStrange::MakeQTPC(AliESDEvent *tESD) {
2317 if(fUseFP) fTPCevent->ClearFast(); // flow package
2318 fQTPCACos=fQTPCASin=fQTPCA=0;
2319 fQTPCCCos=fQTPCCSin=fQTPCC=0;
2323 Double_t eta, phi, w;
2326 Int_t rawN = tESD->GetNumberOfTracks();
2327 for(Int_t id=0; id!=rawN; ++id) {
2328 track = tESD->GetTrack(id);
2330 if( fExcludeTPCEdges )
2331 if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tESD->GetMagneticField() ) ) continue;
2332 if(!PassesFilterBit(track)) continue;
2333 if(!PassesRFPTPCCuts(track)) continue;
2339 fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2340 fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2342 fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
2344 fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2345 fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2347 fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
2350 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2351 if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flowpackage
2354 fTPCevent->SetNumberOfRPs(rfp);
2355 if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
2358 //=======================================================================
2359 void AliAnalysisTaskFlowStrange::AddMCParticleSpy(TList *me) {
2363 tH1D = new TH1D("Pt", "Pt", fPtBins,fPtBinEdge); me->Add(tH1D);
2364 tH1D = new TH1D("Phi", "Phi", 100,0,TMath::TwoPi()); me->Add(tH1D);
2365 tH1D = new TH1D("Eta", "Eta", 100,-1,+1); me->Add(tH1D);
2366 tH1D = new TH1D("Y", "Y", 100,-1,+1); me->Add(tH1D);
2367 tH1D = new TH1D("Rad2", "Rad2", 1000,0,+100); me->Add(tH1D);
2368 tH2D = new TH2D("Dphi", "phi-MCEP;pt;dphi",fPtBins,fPtBinEdge, 72,0,TMath::Pi()); me->Add(tH2D);
2369 tPro = new TProfile("Cos2dphi","Cos2dphi",fPtBins,fPtBinEdge); me->Add(tPro);
2372 //=======================================================================
2373 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, AliAODMCParticle *p) {
2374 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2375 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2376 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Y" ))->Fill( p->Y() );
2377 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2378 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Xv()*p->Xv() +
2379 p->Yv()*p->Yv() ) );
2380 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Dphi" ))->Fill( p->Pt(), GetMCDPHI(p->Phi()) );
2381 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphi" ))->Fill( p->Pt(), TMath::Cos( 2*GetMCDPHI(p->Phi()) ), 1 );
2384 //=======================================================================
2385 Double_t AliAnalysisTaskFlowStrange::GetMCDPHI(Double_t phi) {
2386 Double_t dDPHI = phi - fMCEP;
2387 //if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
2388 //if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
2391 //=======================================================================
2392 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, TParticle *p) {
2393 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2394 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2395 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2396 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Vx()*p->Vx() +
2397 p->Vy()*p->Vy() ) );
2400 //=======================================================================
2401 void AliAnalysisTaskFlowStrange::AddCandidatesSpy(TList *me,Bool_t res) {
2406 tH2D = new TH2D("PhiEta", "PhiEta;Phi;Eta", 100,0,TMath::TwoPi(),100,-1,+1); me->Add(tH2D);
2407 tH2D = new TH2D("PtRAP", "PtRAP;Pt;Y", fPtBins,fPtBinEdge,100,-1,+1); me->Add(tH2D);
2408 tH2D = new TH2D("PtDCA", "PtDCA;Pt;DCA", fPtBins,fPtBinEdge,100,0,1.5); me->Add(tH2D);
2409 tH2D = new TH2D("PtCTP", "PtCTP;Pt;CTP", fPtBins,fPtBinEdge,100,0.95,+1);me->Add(tH2D);
2410 tH2D = new TH2D("PtD0D0", "PtD0D0;Pt;D0D0", fPtBins,fPtBinEdge,100,-5,+5); me->Add(tH2D);
2411 tH2D = new TH2D("PtRad2", "PtRad2;Pt;RadXY",fPtBins,fPtBinEdge,100,0,+50); me->Add(tH2D);
2412 tH2D = new TH2D("PtDL", "PtDL;Pt;DL", fPtBins,fPtBinEdge,100,0,+50); me->Add(tH2D);
2413 tH2D = new TH2D("PtDLlab", "PtDL;Pt;DLlab", fPtBins,fPtBinEdge,100,0,+100); me->Add(tH2D);
2414 tH2D = new TH2D("PtMASS", "PtMASS;Pt;MASS", fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); me->Add(tH2D);
2415 tH2D = new TH2D("APPOS", "APPOS;alphaPOS;QtPOS",100,-2,+2,100,0,0.3); me->Add(tH2D);
2416 tH2D = new TH2D("D0PD0N", "D0PD0N;D0P;D0N", 200,-10,+10,200,-10,+10); me->Add(tH2D);
2417 tH2D = new TH2D("XPOSXNEG","XPOSXNEG;XPOS;XNEG", 200,-50,+50,200,-50,+50); me->Add(tH2D);
2420 tH1D = new TH1D("MCOrigin","MCOrigin;Rad2", 1000,0,50); me->Add(tH1D);
2421 tH2D = new TH2D("PHIRes","PHIRes;PHI;MC-DAT", 72,0,TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2422 tH2D = new TH2D("ETARes","ETARes;ETA;MC-DAT", 16,-0.8,+0.8,100,-0.2,+0.2); me->Add(tH2D);
2423 tH2D = new TH2D("PTRes", "PTRes;Pt;MC-DAT", fPtBins,fPtBinEdge,100,-0.4,+0.4); me->Add(tH2D);
2424 tH2D = new TH2D("RXYRes","RXYRes;RXY;MC-DAT", 100,0,50,100,-4.0,+4.0); me->Add(tH2D);
2426 tH2D = new TH2D("PTDPHIMC","PtDPHIMC;Pt;PHI-MCEP",fPtBins,fPtBinEdge,72,0,TMath::Pi()); me->Add(tH2D);
2427 tPro = new TProfile("Cos2dphiMC", "Cos2dphiMC",fPtBins,fPtBinEdge); me->Add(tPro);
2428 tPro2=new TProfile2D("C2DPHIMCMASS","C2DPHIMCMASS",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); me->Add(tPro2);
2432 //=======================================================================
2433 void AliAnalysisTaskFlowStrange::FillCandidateSpy(TString listName, Bool_t fillRes) {
2434 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PhiEta"))->Fill( fDecayPhi, fDecayEta );
2435 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRAP" ))->Fill( fDecayPt, fDecayRapidity );
2436 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDCA" ))->Fill( fDecayPt, fDecayDCAdaughters );
2437 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtCTP" ))->Fill( fDecayPt, fDecayCosinePointingAngleXY );
2438 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtD0D0"))->Fill( fDecayPt, fDecayProductIPXY );
2439 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRad2"))->Fill( fDecayPt, fDecayRadXY );
2440 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDL" ))->Fill( fDecayPt, fDecayDecayLength );
2441 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDLlab"))->Fill( fDecayPt, fDecayDecayLengthLab );
2442 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtMASS"))->Fill( fDecayPt, fDecayMass );
2443 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("APPOS" ))->Fill( fDecayAlpha, fDecayQt );
2444 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("D0PD0N"))->Fill( fDecayIPpos, fDecayIPneg );
2445 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("XPOSXNEG"))->Fill( fDecayXpos, fDecayXneg );
2447 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTDPHIMC" ))->Fill( fDecayPt, GetMCDPHI( fDecayPhi ) );
2448 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphiMC" ))->Fill( fDecayPt, TMath::Cos( 2*GetMCDPHI(fDecayPhi) ), 1 );
2449 ((TProfile2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("C2DPHIMCMASS" ))->Fill(fDecayPt,fDecayMass,TMath::Cos(2*GetMCDPHI(fDecayPhi)), 1 );
2451 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("MCOrigin"))->Fill( fDecayMatchOrigin );
2452 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDecayPhi, fDecayMatchPhi-fDecayPhi );
2453 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDecayEta, fDecayMatchEta-fDecayEta );
2454 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDecayPt, fDecayMatchPt-fDecayPt );
2455 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("RXYRes"))->Fill( fDecayRadXY, fDecayMatchRadXY-fDecayRadXY );
2459 //=======================================================================
2460 Bool_t AliAnalysisTaskFlowStrange::AcceptCandidate() {
2461 if(fDecayEta<fDecayMinEta) return kFALSE;
2462 if(fDecayEta>fDecayMaxEta) return kFALSE;
2463 if(fDecayPt<fDecayMinPt) return kFALSE;
2464 if(fDecayProductIPXY>fDecayMaxProductIPXY) return kFALSE;
2465 if(fDecayDCAdaughters>fDecayMaxDCAdaughters) return kFALSE;
2466 if(fDecayCosinePointingAngleXY<fDecayMinCosinePointingAngleXY) return kFALSE;
2467 if(fDecayRadXY<fDecayMinRadXY) return kFALSE;
2468 if(TMath::Abs(fDecayRapidity)>fDecayMaxRapidity) return kFALSE;
2470 if(fDecayAPCutPie) {
2471 if(fDecayQt/TMath::Abs(fDecayAlpha)<fDecayMinQt) return kFALSE;
2473 if(fDecayQt<fDecayMinQt) return kFALSE;
2475 if(fDecayDecayLength>fDecayMaxDecayLength*2.6842) return kFALSE;
2477 if(fDecayDecayLength>fDecayMaxDecayLength*7.89) return kFALSE;
2479 if(fSpecie==1) if(fDecayAlpha>0) return kFALSE;
2480 if(fSpecie==2) if(fDecayAlpha<0) return kFALSE;
2483 //=======================================================================
2484 void AliAnalysisTaskFlowStrange::AddTrackSpy(TList *me,Bool_t res) {
2486 tH2D = new TH2D("PHIETA", "PHIETA;PHI;ETA", 100,0,TMath::TwoPi(),100,-2,2); me->Add(tH2D);
2487 tH2D = new TH2D("PTTRACKDECAY", "PTTRACKDECAY;PT;PT", 100,0,10,fPtBins,fPtBinEdge); me->Add(tH2D);
2488 tH2D = new TH2D("IPXYIPZ", "IPXYIPZ;IPXY;IPZ", 100,-10,+10,100,-10,+10); me->Add(tH2D);
2489 tH2D = new TH2D("PTTPCNCLS", "PTTPCNCLS;PT;NCLS", fPtBins,fPtBinEdge,170,0,170); me->Add(tH2D);
2490 tH2D = new TH2D("PTITSNCLS", "PTITSNCLS;PT;NCLS", fPtBins,fPtBinEdge,7,-0.5,6.5); me->Add(tH2D);
2491 tH2D = new TH2D("PTITSLAY", "PTITSLAY;PT;ITSLAYER", fPtBins,fPtBinEdge,6,-0.5,+5.5);me->Add(tH2D);
2492 tH2D = new TH2D("PTITSTPCrefit","PTITSTPCrefit;PT", fPtBins,fPtBinEdge,2,-0.5,+1.5);me->Add(tH2D);
2493 tH2D->GetYaxis()->SetBinLabel(1,"ITS refit");
2494 tH2D->GetYaxis()->SetBinLabel(2,"TPC refit");
2495 tH2D = new TH2D("POSTPCNCLCHI2","POSTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8); me->Add(tH2D);
2496 tH2D = new TH2D("POSTPCNFCLNXR","POSTPCNFCLNXR;NFCLS;NXR", 170,0,170,170,0,170); me->Add(tH2D);
2497 tH2D = new TH2D("POSTPCNCLNFCL","POSTPCNCLNFCL;NCLS;NFCLS", 170,0,170,170,0,170); me->Add(tH2D);
2498 tH2D = new TH2D("POSTPCNCLNSCL","POSTPCNCLNSCL;NCLS;NSCLS", 170,0,170,170,0,170); me->Add(tH2D);
2499 tH2D = new TH2D("NEGTPCNCLCHI2","NEGTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8); me->Add(tH2D);
2500 tH2D = new TH2D("NEGTPCNFCLNXR","NEGTPCNFCLNXR;NFCLS;NXR", 170,0,170,170,0,170); me->Add(tH2D);
2501 tH2D = new TH2D("NEGTPCNCLNFCL","NEGTPCNCLNFCL;NCLS;NFCLS", 170,0,170,170,0,170); me->Add(tH2D);
2502 tH2D = new TH2D("NEGTPCNCLNSCL","NEGTPCNCLNSCL;NCLS;NSCLS", 170,0,170,170,0,170); me->Add(tH2D);
2505 tPro = new TProfile("COSNDPHIMC","COSNDPHIMC",fPtBins,fPtBinEdge); me->Add(tPro);
2508 tH2D = new TH2D("PHIRes", "PHIRes;PHI;MC-DAT", 72,0,TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2509 tH2D = new TH2D("ETARes", "ETARes;ETA;MC-DAT", 16,-0.8,+0.8,100,-0.2,+0.2); me->Add(tH2D);
2510 tH2D = new TH2D("PTRes", "PTRes;Pt;MC-DAT", fPtBins,fPtBinEdge,100,-0.4,+0.4); me->Add(tH2D);
2513 //=======================================================================
2514 void AliAnalysisTaskFlowStrange::FillTrackSpy(TString listName,Bool_t res) {
2515 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PHIETA" ))->Fill( fDaughterPhi, fDaughterEta );
2516 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTRACKDECAY" ))->Fill( fDaughterPt, fDecayPt );
2517 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "IPXYIPZ" ))->Fill( fDaughterImpactParameterXY, fDaughterImpactParameterZ );
2518 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTPCNCLS" ))->Fill( fDaughterPt, fDaughterNClsTPC );
2519 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSNCLS" ))->Fill( fDaughterPt, fDaughterNClsITS );
2520 if( TESTBIT(fDaughterITScm,0) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 0 );
2521 if( TESTBIT(fDaughterITScm,1) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 1 );
2522 if( TESTBIT(fDaughterITScm,2) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 2 );
2523 if( TESTBIT(fDaughterITScm,3) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 3 );
2524 if( TESTBIT(fDaughterITScm,4) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 4 );
2525 if( TESTBIT(fDaughterITScm,5) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 5 );
2526 if( (fDaughterStatus&AliESDtrack::kITSrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 0 );
2527 if( (fDaughterStatus&AliESDtrack::kTPCrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 1 );
2529 if(fDaughterCharge>0) ch="POS";
2530 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLCHI2",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterChi2PerNClsTPC );
2531 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNFCLNXR",ch.Data()) ))->Fill( fDaughterNFClsTPC, fDaughterXRows );
2532 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNFCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNFClsTPC );
2533 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNSCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNSClsTPC );
2535 Double_t cosn = TMath::Cos( fHarmonic*GetMCDPHI(fDaughterPhi) );
2536 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("COSNDPHIMC" ))->Fill( fDaughterPt, cosn, 1 );
2539 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDaughterPhi, fDaughterMatchPhi-fDaughterAtSecPhi );
2540 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDaughterEta, fDaughterMatchEta-fDaughterAtSecEta );
2541 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDaughterPt, fDaughterMatchPt-fDaughterAtSecPt );
2544 //=======================================================================
2545 void AliAnalysisTaskFlowStrange::LoadTrack(AliESDtrack *myTrack, Double_t aodChi2NDF) {
2546 fDaughterCharge = myTrack->Charge();
2547 fDaughterXRows = myTrack->GetTPCCrossedRows();
2548 fDaughterNFClsTPC = myTrack->GetTPCNclsF();
2549 fDaughterNSClsTPC = myTrack->GetTPCnclsS();
2550 fDaughterNClsTPC = myTrack->GetTPCclusters(0);
2552 if(fDaughterNClsTPC>0) fDaughterChi2PerNClsTPC = myTrack->GetTPCchi2()/fDaughterNClsTPC;
2554 fDaughterChi2PerNClsTPC = aodChi2NDF;
2556 myTrack->GetImpactParameters(fDaughterImpactParameterXY,fDaughterImpactParameterZ);
2557 fDaughterStatus = myTrack->GetStatus();
2558 fDaughterITScm = myTrack->GetITSClusterMap();
2559 fDaughterPhi = myTrack->Phi();
2560 fDaughterEta = myTrack->Eta();
2561 fDaughterPt = myTrack->Pt();
2562 fDaughterKinkIndex = myTrack->GetKinkIndex(0);
2564 for(Int_t lay=0; lay!=6; ++lay)
2565 if(TESTBIT(fDaughterITScm,lay)) fDaughterNClsITS++;
2567 //=======================================================================
2568 Bool_t AliAnalysisTaskFlowStrange::AcceptDaughter(Bool_t strongITS) {
2569 if(fDaughterKinkIndex>0) return kFALSE;
2570 if( (fDaughterStatus&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2571 if(fDaughterNFClsTPC<1) return kFALSE;
2572 if(fDaughterPt<fDaughterMinPt) return kFALSE;
2573 if(fDaughterEta<fDaughterMinEta) return kFALSE;
2574 if(fDaughterEta>fDaughterMaxEta) return kFALSE;
2575 if(fDaughterNClsTPC<fDaughterMinNClsTPC) return kFALSE;
2576 if(fDaughterXRows<fDaughterMinXRows) return kFALSE;
2577 if(fDaughterChi2PerNClsTPC>fDaughterMaxChi2PerNClsTPC) return kFALSE;
2578 if(TMath::Abs(fDaughterImpactParameterXY)<fDaughterMinImpactParameterXY) return kFALSE;
2579 if(fDaughterXRows<fDaughterMinXRowsOverNClsFTPC*fDaughterNFClsTPC) return kFALSE;
2581 if( (fDaughterITSrefit) & ((fDaughterStatus&AliESDtrack::kITSrefit)==0) ) return kFALSE;
2582 for(Int_t lay=0; lay!=6; ++lay)
2583 if(fDaughterITSConfig[lay]>-0.5) {
2584 if(fDaughterITSConfig[lay]) {
2585 if(!TESTBIT(fDaughterITScm,lay)) return kFALSE;
2587 if(TESTBIT(fDaughterITScm,lay)) return kFALSE;
2590 if(fDaughterNClsITS<fDaughterMinNClsITS) return kFALSE;
2591 if(fDaughterSPDRequireAny) {
2592 if( !TESTBIT(fDaughterITScm,0)&&!TESTBIT(fDaughterITScm,1)) return kFALSE;
2597 //=======================================================================
2598 Double_t AliAnalysisTaskFlowStrange::GetWDist(const AliVVertex* v0, const AliVVertex* v1) {
2599 // calculate sqrt of weighted distance to other vertex
2601 printf("One of vertices is not valid\n");
2604 static TMatrixDSym vVb(3);
2606 double dx = v0->GetX()-v1->GetX();
2607 double dy = v0->GetY()-v1->GetY();
2608 double dz = v0->GetZ()-v1->GetZ();
2609 double cov0[6],cov1[6];
2610 v0->GetCovarianceMatrix(cov0);
2611 v1->GetCovarianceMatrix(cov1);
2612 vVb(0,0) = cov0[0]+cov1[0];
2613 vVb(1,1) = cov0[2]+cov1[2];
2614 vVb(2,2) = cov0[5]+cov1[5];
2615 vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2616 vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2618 if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2619 dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2620 + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2621 return dist>0 ? TMath::Sqrt(dist) : -1;
2623 //=======================================================================
2624 Bool_t AliAnalysisTaskFlowStrange::plpMV(const AliVEvent *event) {
2625 // check for multi-vertexer pile-up
2626 const AliAODEvent *aod = (const AliAODEvent*)event;
2627 const AliESDEvent *esd = (const AliESDEvent*)event;
2629 const int kMinPlpContrib = 5;
2630 const double kMaxPlpChi2 = 5.0;
2631 const double kMinWDist = 15;
2634 printf("Event is neither of AOD nor ESD\n");
2638 const AliVVertex* vtPrm = 0;
2639 const AliVVertex* vtPlp = 0;
2643 if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2644 vtPrm = aod->GetPrimaryVertex();
2645 if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2648 if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
2649 vtPrm = esd->GetPrimaryVertexTracks();
2650 if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
2652 //int bcPrim = vtPrm->GetBC();
2654 for (int ipl=0;ipl<nPlp;ipl++) {
2655 vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
2657 if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2658 if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2659 // int bcPlp = vtPlp->GetBC();
2660 // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2662 double wDst = GetWDist(vtPrm,vtPlp);
2663 if (wDst<kMinWDist) continue;
2665 return kTRUE; // pile-up: well separated vertices
2671 //=======================================================================
2672 void AliAnalysisTaskFlowStrange::MakeFilterBits() {
2674 fFB1 = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
2676 fFB1024 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
2677 fFB1024->SetMinNCrossedRowsTPC(120);
2678 fFB1024->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
2679 fFB1024->SetMaxChi2PerClusterITS(36);
2680 fFB1024->SetMaxFractionSharedTPCClusters(0.4);
2681 fFB1024->SetMaxChi2TPCConstrainedGlobal(36);
2682 fFB1024->SetEtaRange(-0.9,0.9);
2683 fFB1024->SetPtRange(0.15, 1e10);
2685 //=======================================================================
2686 Bool_t AliAnalysisTaskFlowStrange::PassesFilterBit(AliESDtrack *track) {
2688 switch(fRFPFilterBit) {
2690 ret = fFB1024->AcceptTrack(track);
2693 ret = fFB1->AcceptTrack(track);
2697 //=======================================================================
2698 void AliAnalysisTaskFlowStrange::LoadVZEROResponse() {
2700 TString run = fVZEResponse->GetTitle();
2701 if( run.Atoi() == fRunNumber ) return;
2702 fVZEResponse = NULL;
2705 fVZEResponse = dynamic_cast<TH2D*> (fVZEload->FindObject( Form("%d",fRunNumber) ));
2707 printf("New VZE calibration: run %d || %s -> Entries %.0f\n",fRunNumber, fVZEResponse->GetTitle(),fVZEResponse->GetEntries());
2708 //=>external weights
2709 for(int i=0;i!=64;++i) fVZEextW[i]=1;
2711 Double_t minC = fCentPerMin, maxC = fCentPerMax;
2716 Int_t ybinmin = fVZEResponse->GetYaxis()->FindBin(minC+1e-6);
2717 Int_t ybinmax = fVZEResponse->GetYaxis()->FindBin(maxC-1e-6);
2718 if(fSkipCentralitySelection) {
2722 for(int i=0;i!=64;++i) fVZEextW[i] = fVZEResponse->Integral(i+1,i+1,ybinmin,ybinmax)/(maxC-minC);
2723 //ring-wise normalization
2725 for(int j=0; j!=8; ++j) {
2727 for(int i=0;i!=8;++i) ring[j] += fVZEextW[j*8+i]/8;
2729 //disk-wise normalization
2731 int xbinmin, xbinmax;
2732 xbinmin = 1+8*fVZECa;
2733 xbinmax = 8+8*fVZECb;
2734 disk[0] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
2735 xbinmin = 33+8*fVZEAa;
2736 xbinmax = 40+8*fVZEAb;
2737 disk[1] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
2738 //for(int i=0;i!=64;++i) printf("CELL %d -> W = %f ||",i,fVZEextW[i]);
2740 for(int i=0;i!=64;++i) fVZEextW[i] = disk[i/32]/fVZEextW[i];
2742 for(int i=0;i!=64;++i) fVZEextW[i] = ring[i/8]/fVZEextW[i];
2744 //for(int i=0;i!=64;++i) printf(" W = %f \n",fVZEextW[i]);
2747 printf("VZE calibration: requested but not found!!!\n");
2750 //=======================================================================
2751 void AliAnalysisTaskFlowStrange::AddVZEQA() {
2754 TList *tList = new TList();
2755 tList->SetName( "VZEQA" );
2757 fList->Add( tList );
2758 tH2D = new TH2D("EQU","EQU;VZEeqmult-VZEmult;cell",100,-5,+5,64,0,64); tList->Add( tH2D );
2759 prof = new TProfile2D("LINbefCAL","LINbef;VZEcell;VZEeqmult;SPDtrkl", 64,0,64,350,0,700,0,10000); tList->Add( prof );
2760 prof = new TProfile2D("LINaftCAL","LINaft;VZEcell;VZEeqmult;SPDtrkl", 64,0,64,350,0,700,0,10000); tList->Add( prof );
2762 //=======================================================================
2763 void AliAnalysisTaskFlowStrange::FillVZEQA() {
2764 AliESDEvent *tESD = (AliESDEvent*) (InputEvent());
2765 AliAODEvent *tAOD = (AliAODEvent*) (InputEvent());
2774 //=======================================================================
2775 void AliAnalysisTaskFlowStrange::FillVZEQA(AliAODEvent *tAOD) {
2776 AliVVZERO *vzero = tAOD->GetVZEROData();
2777 AliAODTracklets *tracklets = tAOD->GetTracklets();
2779 if(!tracklets) return;
2780 Double_t mult, eqmult;
2781 Int_t trkl = tracklets->GetNumberOfTracklets();
2782 for(int id=0; id!=64; ++id) {
2783 mult = vzero->GetMultiplicity(id);
2784 eqmult = tAOD->GetVZEROEqMultiplicity(id);
2785 ((TH2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject("EQU"))->Fill(eqmult-mult,id);
2786 ((TProfile2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject( "LINbefCAL" ))->Fill(id,eqmult,trkl,1);
2787 ((TProfile2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject( "LINaftCAL" ))->Fill(id,eqmult*fVZEextW[id],trkl,1);
2790 //=======================================================================
2791 void AliAnalysisTaskFlowStrange::AddVZEROResponse() {
2792 fVZEResponse = NULL;
2793 AliVEvent *event = InputEvent();
2795 Int_t thisrun = event->GetRunNumber();
2796 fVZEResponse = new TH2D( Form("%d",thisrun), Form("%d;cell;CC",thisrun), 64,0,64, 110, -10, 100);
2797 fList->Add(fVZEResponse);
2799 //=======================================================================
2800 void AliAnalysisTaskFlowStrange::SaveVZEROResponse() {
2801 if(!fVZEResponse) return;
2802 AliVEvent *event = InputEvent();
2805 // reject event with low ocupancy in VZERO
2806 // for centralities below 60% this should not happen anyway
2807 Double_t rejectEvent = kFALSE;
2808 for(int id=0; id!=64; ++id) {
2809 w = event->GetVZEROEqMultiplicity(id);
2810 if(w<3) rejectEvent = kTRUE;
2812 if(rejectEvent) return;
2814 for(int id=0; id!=64; ++id) {
2815 w = event->GetVZEROEqMultiplicity(id);
2816 fVZEResponse->Fill(id,fThisCent,w);
2819 //=======================================================================
2820 Int_t AliAnalysisTaskFlowStrange::RefMult(AliAODEvent *tAOD, Int_t fb) {
2822 for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
2823 AliAODTrack *t = dynamic_cast<AliAODTrack*>(tAOD->GetTrack( i ));
2825 if( !t->TestFilterBit(fb) ) continue;
2826 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2827 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2828 if( t->GetTPCNcls()<70 ) continue;
2829 //if( t->GetTPCsignal()<10.0 ) continue;
2830 if( t->Chi2perNDF()<0.2 ) continue;
2835 //=======================================================================
2836 Int_t AliAnalysisTaskFlowStrange::RefMultTPC() {
2837 AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2840 for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2841 AliAODTrack *t = dynamic_cast<AliAODTrack*>(ev->GetTrack( i ));
2843 if( !t->TestFilterBit(1) ) continue;
2844 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2845 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2846 if( t->GetTPCNcls()<70 ) continue;
2847 if( t->GetTPCsignal()<10.0 ) continue;
2848 if( t->Chi2perNDF()<0.2 ) continue;
2853 //=======================================================================
2854 Int_t AliAnalysisTaskFlowStrange::RefMultGlobal() {
2855 AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2858 for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2859 AliAODTrack *t = dynamic_cast<AliAODTrack*>(ev->GetTrack( i ));
2861 if( !t->TestFilterBit(16) ) continue;
2862 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2863 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2864 if( t->GetTPCNcls()<70 ) continue;
2865 if( t->GetTPCsignal()<10.0 ) continue;
2866 if( t->Chi2perNDF()<0.1 ) continue;
2867 Double_t b[3], bcov[3];
2868 if( !t->PropagateToDCA(ev->GetPrimaryVertex(),ev->GetMagneticField(),100,b,bcov) ) continue;
2869 if( b[0]>+0.3 || b[0]<-0.3 || b[1]>+0.3 || b[1]<-0.3) continue;
2874 //=======================================================================
2875 void AliAnalysisTaskFlowStrange::ResetContainers() {
2877 fTPCevent->ClearFast();
2878 fVZEevent->ClearFast();
2881 //=======================================================================
2882 void AliAnalysisTaskFlowStrange::AddTrackVn(TList *tList) {
2886 tProfile = new TProfile("SP_uVZEA","u x Q_{VZEA}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2887 tProfile = new TProfile("SP_uVZEC","u x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2888 tProfile = new TProfile("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2889 tProfile = new TProfile("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2890 tProfile = new TProfile("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2892 tProfile = new TProfile("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2893 tProfile = new TProfile("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2894 tProfile = new TProfile("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2896 tProfile = new TProfile("SP_uTPCA","u x Q_{TPCA}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2897 tProfile = new TProfile("SP_uTPCC","u x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2898 tProfile = new TProfile("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2900 tProfile = new TProfile("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2901 tProfile = new TProfile("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2903 tH1D = new TH1D("QC_HistPt_P","HistPt_P",fPtBins,fPtBinEdge); tList->Add( tH1D );
2904 tH1D = new TH1D("QC_HistPt_Q","HistPt_Q",fPtBins,fPtBinEdge); tList->Add( tH1D );
2906 tProfile = new TProfile("QC_C2","QC_C2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2907 tProfile = new TProfile("QC_C4","QC_C4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2908 tProfile = new TProfile("QC_DC2","QC_DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2909 tProfile = new TProfile("QC_DC4","QC_DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2910 tProfile = new TProfile("QC_C2C4","QC_C2C4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2911 tProfile = new TProfile("QC_C2DC2","QC_C2DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2912 tProfile = new TProfile("QC_C2DC4","QC_C2DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2913 tProfile = new TProfile("QC_C4DC2","QC_C4DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2914 tProfile = new TProfile("QC_C4DC4","QC_C4DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2915 tProfile = new TProfile("QC_DC2DC4","QC_DC2DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2917 tProfile = new TProfile("QC_pCos","QC_pCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2918 tProfile = new TProfile("QC_pSin","QC_pSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2919 tProfile = new TProfile("QC_qCos","QC_qCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2920 tProfile = new TProfile("QC_qSin","QC_qSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2921 tProfile = new TProfile("QC_q2hCos","QC_q2hCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2922 tProfile = new TProfile("QC_q2hSin","QC_q2hSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2924 tH1D = new TH1D("SP_vnVZEA","SP_vnVZEA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2925 tH1D = new TH1D("SP_vnVZEC","SP_vnVZEC",fPtBins,fPtBinEdge); tList->Add( tH1D );
2926 tH1D = new TH1D("SP_vnTPCA","SP_vnTPCA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2927 tH1D = new TH1D("SP_vnTPCC","SP_vnTPCC",fPtBins,fPtBinEdge); tList->Add( tH1D );
2928 tH1D = new TH1D("QC_Cum2","QC_Cum2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2929 tH1D = new TH1D("QC_Cum4","QC_Cum4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2930 tH1D = new TH1D("QC_DCum2","QC_DCum2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2931 tH1D = new TH1D("QC_DCum4","QC_DCum4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2932 tH1D = new TH1D("SP_vnVZEGA","SP_vnVZEGA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2933 tH1D = new TH1D("SP_vnVZEWA","SP_vnVZEWA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2934 tH1D = new TH1D("SP_vnTPCAA","SP_vnTPCAA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2935 tH1D = new TH1D("QC_vn2","QC_vn2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2936 tH1D = new TH1D("QC_vn4","QC_vn4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2939 tProfile = new TProfile("MC_COSNDPHI","MC_COSNDPHI",fPtBins,fPtBinEdge,-3,+3); tList->Add( tProfile );
2940 tH2D = new TH2D("MC_COSNDPHI_uQVZEA","MC_COSNDPHI_uQVZEA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2941 tH2D = new TH2D("MC_COSNDPHI_uQVZEC","MC_COSNDPHI_uQVZEC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2942 tH2D = new TH2D("MC_COSNDPHI_uQTPCA","MC_COSNDPHI_uQTPCA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2943 tH2D = new TH2D("MC_COSNDPHI_uQTPCC","MC_COSNDPHI_uQTPCC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2946 //=======================================================================
2947 void AliAnalysisTaskFlowStrange::AddDecayVn(TList *tList) {
2948 TProfile2D *tProfile;
2951 tH2D = new TH2D("DecayYield_PtMass","Decay_PtMass",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2952 tProfile = new TProfile2D("DecayAvgPt_PtMass","Decay_PtMass",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tProfile ); tProfile->Sumw2();
2954 tProfile = new TProfile2D("SP_uVZEA","u x Q_{VZEA}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2955 tProfile = new TProfile2D("SP_uVZEC","u x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2956 tProfile = new TProfile2D("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2957 tProfile = new TProfile2D("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2958 tProfile = new TProfile2D("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2960 tProfile = new TProfile2D("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2961 tProfile = new TProfile2D("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2962 tProfile = new TProfile2D("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2964 tProfile = new TProfile2D("SP_uTPCA","u x Q_{TPCA}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2965 tProfile = new TProfile2D("SP_uTPCC","u x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2966 tProfile = new TProfile2D("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2968 tProfile = new TProfile2D("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2969 tProfile = new TProfile2D("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2971 tH2D = new TH2D("QC_HistPt_P","HistPt_P",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2972 tH2D = new TH2D("QC_HistPt_Q","HistPt_Q",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2974 tProfile = new TProfile2D("QC_C2","QC_C2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2975 tProfile = new TProfile2D("QC_C4","QC_C4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2976 tProfile = new TProfile2D("QC_DC2","QC_DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2977 tProfile = new TProfile2D("QC_DC4","QC_DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2978 tProfile = new TProfile2D("QC_C2C4","QC_C2C4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2979 tProfile = new TProfile2D("QC_C2DC2","QC_C2DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2980 tProfile = new TProfile2D("QC_C2DC4","QC_C2DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2981 tProfile = new TProfile2D("QC_C4DC2","QC_C4DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2982 tProfile = new TProfile2D("QC_C4DC4","QC_C4DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2983 tProfile = new TProfile2D("QC_DC2DC4","QC_DC2DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2985 tProfile = new TProfile2D("QC_pCos","QC_pCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2986 tProfile = new TProfile2D("QC_pSin","QC_pSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2987 tProfile = new TProfile2D("QC_qCos","QC_qCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2988 tProfile = new TProfile2D("QC_qSin","QC_qSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2989 tProfile = new TProfile2D("QC_q2hCos","QC_q2hCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2990 tProfile = new TProfile2D("QC_q2hSin","QC_q2hSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2992 tH2D = new TH2D("SP_vnVZEA","SP_vnVZEA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2993 tH2D = new TH2D("SP_vnVZEC","SP_vnVZEC",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2994 tH2D = new TH2D("SP_vnTPCA","SP_vnTPCA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2995 tH2D = new TH2D("SP_vnTPCC","SP_vnTPCC",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2996 tH2D = new TH2D("QC_Cum2","QC_Cum2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2997 tH2D = new TH2D("QC_Cum4","QC_Cum4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2998 tH2D = new TH2D("QC_DCum2","QC_DCum2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2999 tH2D = new TH2D("QC_DCum4","QC_DCum4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3000 tH2D = new TH2D("SP_vnVZEGA","SP_vnVZEGA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3001 tH2D = new TH2D("SP_vnVZEWA","SP_vnVZEWA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3002 tH2D = new TH2D("SP_vnTPCAA","SP_vnTPCAA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3003 tH2D = new TH2D("QC_vn2","QC_vn2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3004 tH2D = new TH2D("QC_vn4","QC_vn4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3006 tProfile = new TProfile2D("MC_COSNDPHI","MC_COSNDPHI",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tProfile );
3009 //=======================================================================
3010 TList* AliAnalysisTaskFlowStrange::RebinDecayVn(Int_t nbins, Int_t *bins) {
3011 TList *out = new TList();
3013 out->SetName( fList->GetName() );
3017 ori = (TList*) fList->FindObject("Event");
3020 ori = (TList*) fList->FindObject("MakeQSpy");
3023 ori = (TList*) fList->FindObject("V0SAllVn");
3024 end = RebinDecayVn(ori,nbins,bins);
3025 end->SetName( ori->GetName() ); out->Add( end );
3027 ori = (TList*) fList->FindObject("V0SSelVn");
3028 end = RebinDecayVn(ori,nbins,bins);
3029 end->SetName( ori->GetName() ); out->Add( end );
3032 ori = (TList*) fList->FindObject("V0SMthVn");
3033 end = RebinDecayVn(ori,nbins,bins);
3034 end->SetName( ori->GetName() ); out->Add( end );
3036 ori = (TList*) fList->FindObject("V0SMthPosPosVn");
3037 end = RebinDecayVn(ori,nbins,bins);
3038 end->SetName( ori->GetName() ); out->Add( end );
3040 ori = (TList*) fList->FindObject("V0SMthNegNegVn");
3041 end = RebinDecayVn(ori,nbins,bins);
3042 end->SetName( ori->GetName() ); out->Add( end );
3044 ori = (TList*) fList->FindObject("V0SMthPosNegVn");
3045 end = RebinDecayVn(ori,nbins,bins);
3046 end->SetName( ori->GetName() ); out->Add( end );
3048 ori = (TList*) fList->FindObject("V0SUnMthVn");
3049 end = RebinDecayVn(ori,nbins,bins);
3050 end->SetName( ori->GetName() ); out->Add( end );
3055 //=======================================================================
3056 TList* AliAnalysisTaskFlowStrange::RebinDecayVn(TList *tList,Int_t nbins, Int_t *bins) {
3057 // getting expected number of mass bins
3059 for(int i=0; i!=nbins; ++i) sum += bins[i];
3061 TList *list = new TList();
3065 Double_t pt[200], mass[200];
3068 tH2D = ((TH2D*)tList->FindObject( "DecayYield_PtMass" ));
3069 list->Add(tH2D); //keeping it as it is
3070 int nmassbins = tH2D->GetNbinsY();
3072 if( nmassbins!=sum ) {
3073 printf("Error: incompatible binning %d vs %d\nBYE\n",nmassbins,sum);
3077 npt = tH2D->GetNbinsX();
3078 for(int i=0; i!=npt+1; ++i) pt[i] = tH2D->GetXaxis()->GetBinLowEdge(i+1);
3080 for(int i=0,j=0; i!=nbins+1; ++i) {
3081 mass[i] = tH2D->GetYaxis()->GetBinLowEdge(j+1);
3084 //TProfile2D migrating info
3085 TProfile2D *tProfileOld, *tProfileNew;
3086 TString tprofiles[31] = {"DecayAvgPt_PtMass","SP_uVZEA","SP_uVZEC","SP_VZEAVZEC","SP_VZEATPC",
3087 "SP_VZECTPC","SP_uVZEAuVZEC","SP_uVZEAVZEAVZEC","SP_uVZECVZEAVZEC","SP_uTPCA",
3088 "SP_uTPCC","SP_TPCATPCC","SP_uTPCATPCATPCC","SP_uTPCCTPCATPCC","QC_C2",
3089 "QC_C4","QC_DC2","QC_DC4","QC_C2C4","QC_C2DC2",
3090 "QC_C2DC4","QC_C4DC2","QC_C4DC4","QC_DC2DC4","QC_pCos",
3091 "QC_pSin","QC_qCos","QC_qSin","QC_q2hCos","QC_q2hSin",
3093 for(int i=0; i!=31; ++i) {
3094 tProfileOld = (TProfile2D*) tList->FindObject( tprofiles[i].Data() );
3095 if(!tProfileOld) continue;
3096 TArrayD *oldsw2 = tProfileOld->GetSumw2();
3097 TArrayD *oldbsw2 = tProfileOld->GetBinSumw2();
3098 tProfileNew = new TProfile2D( tprofiles[i].Data(), tprofiles[i].Data(),
3099 npt,pt,nbins,mass,"s");
3100 tProfileNew->Sumw2();
3101 list->Add( tProfileNew );
3102 TArrayD *newsw2 = tProfileNew->GetSumw2();
3103 TArrayD *newbsw2 = tProfileNew->GetBinSumw2();
3104 for(int x=1; x!=tProfileOld->GetNbinsX()+1; ++x) { //pt
3105 for(int y=1; y!=tProfileNew->GetNbinsY()+1; ++y) { //mass
3106 Double_t sContent=0;
3107 Double_t sEntries=0;
3108 Double_t sSqWeigh=0;
3109 Double_t sSqBWeig=0;
3110 Int_t binnew = tProfileNew->GetBin(x,y);
3111 Double_t minmass = tProfileNew->GetYaxis()->GetBinLowEdge( y );
3112 Double_t maxmass = tProfileNew->GetYaxis()->GetBinLowEdge( y+1 );
3113 Int_t minbin = tProfileOld->GetYaxis()->FindBin( minmass+1e-10 );
3114 Int_t maxbin = tProfileOld->GetYaxis()->FindBin( maxmass-1e-10 );
3115 for(int k=minbin; k!=maxbin+1; ++k) {
3116 Int_t binold = tProfileOld->GetBin(x,k);
3117 Double_t wk = tProfileOld->GetBinEntries( binold );
3118 Double_t yk = tProfileOld->GetBinContent( binold );
3121 sSqWeigh += oldsw2->At(binold);
3122 if(oldbsw2->GetSize()) sSqBWeig += oldbsw2->At(binold);
3124 tProfileNew->SetBinEntries( binnew, sEntries );
3125 tProfileNew->SetBinContent( binnew, sContent );
3126 newsw2->SetAt(sSqWeigh, binnew);
3127 if(oldbsw2->GetSize()) newbsw2->SetAt(sSqBWeig, binnew);
3132 //TH2D all new blank
3133 TString th2ds[15] = {"QC_HisPt_P","QC_HistPt_Q","SP_vnVZEA","SP_vnVZEC","SP_vnTPCA",
3134 "SP_vnTPCC","QC_Cum2","QC_Cum4","QC_DCum2","QC_DCum4",
3135 "SP_vnVZEGA","SP_vnVZEWA","SP_vnTPCAA","QC_vn2","QC_vn4"};
3136 for(int i=0; i!=15; ++i) {
3137 tH2D = new TH2D( th2ds[i].Data(),th2ds[i].Data(),
3143 //int nmass = tH2Dold->GetNbinsY();
3148 for(int i=0; i!=38; ++i) {
3149 tProfile = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( nam[i].Data() ));
3150 if(i==o) { // binning
3151 nxs = tProfile->GetXaxis()->GetNbins();
3152 nys = tProfile->GetYaxis()->GetNbins();
3156 for(int y=0; y!=nys; ++y)
3157 ny[y] = tProfile->GetYaxis()->GetBinLowEdge(y+1);
3161 for(int x=0; x!=nxs; ++x)
3162 nx[x] = tProfile->GetXaxis()->GetBinLowEdge(x+1);
3165 tProfileNew = new TProfile2D(tProfile->GetName(),tProfile->GetTitle(),nxs,nx,nys,ny,"s"); list->Add( tProfileNew );
3167 for(int y=0; y!=nys; ++y) {
3169 for(int x=0; x!=nxs; ++x) {
3178 //=======================================================================
3179 void AliAnalysisTaskFlowStrange::QCStoreTrackVn(TString name) {
3180 // getting transients
3181 TProfile *pCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ));
3182 TProfile *pSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ));
3183 TProfile *qCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ));
3184 TProfile *qSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ));
3185 TProfile *q2hCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ));
3186 TProfile *q2hSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ));
3188 printf("**QCStoreTrackVn( %s )\n",name.Data());
3189 printf(" Read %.0f entries in p set and %.0f entries in q set\n",pCos->GetEntries(),qCos->GetEntries());
3190 printf(" pCos[5] %.16f | pSin[5] %.16f \n", pCos->GetBinContent(5), pSin->GetBinContent(5));
3191 printf(" qCos[5] %.16f | qSin[5] %.16f \n", qCos->GetBinContent(5), qSin->GetBinContent(5));
3192 printf(" q2hCos[5] %.16f | q2hSin[5] %.16f \n", q2hCos->GetBinContent(5), q2hSin->GetBinContent(5));
3194 // filling {2} and {4} correlator
3195 TProfile *c2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
3196 TProfile *c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
3197 TProfile *dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
3198 TProfile *dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
3199 TProfile *c2c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
3200 TProfile *c2dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
3201 TProfile *c2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
3202 TProfile *c4dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
3203 TProfile *c4dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
3204 TProfile *dc2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
3205 double tpc_qcos = fQTPCACos + fQTPCCCos;
3206 double tpc_qsin = fQTPCASin + fQTPCCSin;
3207 double tpc_qsqr = tpc_qcos*tpc_qcos + tpc_qsin*tpc_qsin;
3208 double tpc_q2hsqr = fQTPC2hCos*fQTPC2hCos + fQTPC2hSin*fQTPC2hSin;
3209 double tpc_qmul = fQTPCA + fQTPCC;
3210 int n = c2->GetNbinsX();
3211 for(int i=1; i!=n+1; ++i) {
3212 double mp = pCos->GetBinEntries( i );
3213 if( mp<1 ) { if(fDebug>2) printf(" bin %d:: mp (%.16f) < 1!\n",i,mp); continue; }
3214 double mm1 = tpc_qmul*(tpc_qmul-1);
3215 if( mm1<1e-100 ) { if(fDebug>2) printf(" bin %d:: mm1<1e-100!\n",i); continue; }
3216 double mq = qCos->GetBinEntries( i );
3217 double mpmmq = mp*tpc_qmul-mq;
3218 if( mpmmq<1e-100 ) { if(fDebug>2) printf(" bin %d:: mpmmq<1e-100!\n",i); continue; }
3219 double pt = c2->GetBinCenter( i );
3220 double pcos = pCos->GetBinContent( i )*mp;
3221 double psin = pSin->GetBinContent( i )*mp;
3222 double qcos = qCos->GetBinContent( i )*mq;
3223 double qsin = qSin->GetBinContent( i )*mq;
3224 double q2hcos = q2hCos->GetBinContent( i )*mq;
3225 double q2hsin = q2hSin->GetBinContent( i )*mq;
3226 double pQ = pcos*tpc_qcos+psin*tpc_qsin;
3227 double q2nQnQn = (qcos*tpc_qcos + qsin*tpc_qsin)*tpc_qcos + (qsin*tpc_qcos-qcos*tpc_qsin)*tpc_qsin;
3228 double pnQnQ2n = (pcos*tpc_qcos - psin*tpc_qsin)*fQTPC2hCos + (psin*tpc_qcos+pcos*tpc_qsin)*fQTPC2hSin;
3229 double tC2 = (tpc_qsqr-tpc_qmul)/mm1;
3230 double tDC2 = (pQ-mq)/mpmmq;
3231 c2->Fill( pt, tC2, mm1 );
3232 dc2->Fill( pt, tDC2, mpmmq );
3233 c2dc2->Fill( pt, tC2*tDC2, mm1*mpmmq );
3234 double mm1m2m3 = tpc_qmul*(tpc_qmul-1)*(tpc_qmul-2)*(tpc_qmul-3);
3235 if(mm1m2m3<1e-100) continue;
3236 double mpm3mqm1m2 = (mp*tpc_qmul-3*mq)*(tpc_qmul-1)*(tpc_qmul-2);
3237 if(mpm3mqm1m2<1e-100) continue;
3238 double req2hqnqn = fQTPC2hCos*(tpc_qcos*tpc_qcos-tpc_qsin*tpc_qsin)+2*fQTPC2hSin*tpc_qcos*tpc_qsin;
3239 double tC4 = (tpc_qsqr*tpc_qsqr + tpc_q2hsqr - 2*req2hqnqn - 2*(2*tpc_qsqr*(tpc_qmul-2)-tpc_qmul*(tpc_qmul-3)))/mm1m2m3;
3240 double tDC4 = pQ*tpc_qsqr -q2nQnQn -pnQnQ2n -2*(tpc_qmul-1)*pQ -2*mq*(tpc_qsqr-tpc_qmul+3) +6*(qcos*tpc_qcos+qsin*tpc_qsin) +(q2hcos*fQTPC2hCos+q2hsin*fQTPC2hSin);
3242 c4->Fill( pt, tC4, mm1m2m3 );
3243 dc4->Fill( pt, tDC4, mpm3mqm1m2 );
3244 c2c4->Fill( pt, tC2*tC4, mm1*mm1m2m3 );
3245 c2dc4->Fill( pt, tC2*tDC4, mm1*mpm3mqm1m2 );
3246 c4dc2->Fill( pt, tC4*tDC2, mm1m2m3*mpmmq );
3247 c4dc4->Fill( pt, tC4*tDC4, mm1m2m3*mpm3mqm1m2 );
3248 dc2dc4->Fill( pt, tDC2*tDC4, mpmmq*mpm3mqm1m2 );
3250 // clean for next event
3251 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Reset();
3252 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Reset();
3253 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Reset();
3254 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Reset();
3255 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Reset();
3256 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Reset();
3258 //=======================================================================
3259 void AliAnalysisTaskFlowStrange::QCStoreDecayVn(TString name) {
3260 // getting transients
3261 TProfile2D *pCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ));
3262 TProfile2D *pSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ));
3263 TProfile2D *qCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ));
3264 TProfile2D *qSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ));
3265 TProfile2D *q2hCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ));
3266 TProfile2D *q2hSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ));
3268 printf("**QCStoreTrackVn( %s )\n",name.Data());
3269 printf(" Read %.0f entries in p set and %.0f entries in q set\n",pCos->GetEntries(),qCos->GetEntries());
3270 printf(" pCos[5][5] %.16f | pSin[5][5] %.16f \n", pCos->GetBinContent(5,5), pSin->GetBinContent(5,5));
3271 printf(" qCos[5][5] %.16f | qSin[5][5] %.16f \n", qCos->GetBinContent(5,5), qSin->GetBinContent(5,5));
3272 printf(" q2hCos[5][5] %.16f | q2hSin[5][5] %.16f \n", q2hCos->GetBinContent(5,5), q2hSin->GetBinContent(5,5));
3274 // filling {2} and {4} correlator
3275 TProfile2D *c2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
3276 TProfile2D *c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
3277 TProfile2D *dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
3278 TProfile2D *dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
3279 TProfile2D *c2c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
3280 TProfile2D *c2dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
3281 TProfile2D *c2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
3282 TProfile2D *c4dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
3283 TProfile2D *c4dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
3284 TProfile2D *dc2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
3285 double tpc_qcos = fQTPCACos + fQTPCCCos;
3286 double tpc_qsin = fQTPCASin + fQTPCCSin;
3287 double tpc_qsqr = tpc_qcos*tpc_qcos + tpc_qsin*tpc_qsin;
3288 double tpc_q2hsqr = fQTPC2hCos*fQTPC2hCos + fQTPC2hSin*fQTPC2hSin;
3289 double tpc_qmul = fQTPCA + fQTPCC;
3290 int n = c2->GetNbinsX();
3291 int m = c2->GetNbinsY();
3292 for(int i=1; i!=n+1; ++i) {
3293 double pt = c2->GetXaxis()->GetBinCenter( i );
3294 for(int j=1; j!=m+1; ++j) {
3295 double ms = c2->GetYaxis()->GetBinCenter( j );
3296 int k = pCos->GetBin(i,j);
3297 double mp = pCos->GetBinEntries( k );
3298 if( mp<1 ) { if(fDebug>2) printf(" bin %d,%d:: mp (%.16f) < 1!\n",i,j,mp); continue; }
3299 double mm1 = tpc_qmul*(tpc_qmul-1);
3300 if( mm1<1e-100 ) { if(fDebug>2) printf(" bin %d,%d:: mm1<1e-100!\n",i,j); continue; }
3301 double mq = qCos->GetBinEntries( k );
3302 double mpmmq = mp*tpc_qmul-mq;
3303 if( mpmmq<1e-100 ) { if(fDebug>2) printf(" bin %d,%d:: mpmmq<1e-100!\n",i,j); continue; }
3304 double pcos = pCos->GetBinContent( i,j )*mp;
3305 double psin = pSin->GetBinContent( i,j )*mp;
3306 double qcos = qCos->GetBinContent( i,j )*mq;
3307 double qsin = qSin->GetBinContent( i,j )*mq;
3308 double q2hcos = q2hCos->GetBinContent( i,j )*mq;
3309 double q2hsin = q2hSin->GetBinContent( i,j )*mq;
3310 double pQ = pcos*tpc_qcos+psin*tpc_qsin;
3311 double q2nQnQn = (qcos*tpc_qcos + qsin*tpc_qsin)*tpc_qcos + (qsin*tpc_qcos-qcos*tpc_qsin)*tpc_qsin;
3312 double pnQnQ2n = (pcos*tpc_qcos - psin*tpc_qsin)*fQTPC2hCos + (psin*tpc_qcos+pcos*tpc_qsin)*fQTPC2hSin;
3313 double tC2 = (tpc_qsqr-tpc_qmul)/mm1;
3314 double tDC2 = (pQ-mq)/mpmmq;
3315 c2->Fill( pt, ms, tC2, mm1 );
3316 dc2->Fill( pt, ms, tDC2, mpmmq );
3317 c2dc2->Fill( pt, ms, tC2*tDC2, mm1*mpmmq );
3318 double mm1m2m3 = tpc_qmul*(tpc_qmul-1)*(tpc_qmul-2)*(tpc_qmul-3);
3319 if(mm1m2m3<1e-100) continue;
3320 double mpm3mqm1m2 = (mp*tpc_qmul-3*mq)*(tpc_qmul-1)*(tpc_qmul-2);
3321 if(mpm3mqm1m2<1e-100) continue;
3322 double req2hqnqn = fQTPC2hCos*(tpc_qcos*tpc_qcos-tpc_qsin*tpc_qsin)+2*fQTPC2hSin*tpc_qcos*tpc_qsin;
3323 double tC4 = (tpc_qsqr*tpc_qsqr + tpc_q2hsqr - 2*req2hqnqn - 2*(2*tpc_qsqr*(tpc_qmul-2)-tpc_qmul*(tpc_qmul-3)))/mm1m2m3;
3324 double tDC4 = pQ*tpc_qsqr -q2nQnQn -pnQnQ2n -2*(tpc_qmul-1)*pQ -2*mq*(tpc_qsqr-tpc_qmul+3) +6*(qcos*tpc_qcos+qsin*tpc_qsin) +(q2hcos*fQTPC2hCos+q2hsin*fQTPC2hSin);
3326 c4->Fill( pt, ms, tC4, mm1m2m3 );
3327 dc4->Fill( pt, ms, tDC4, mpm3mqm1m2 );
3328 c2c4->Fill( pt, ms, tC2*tC4, mm1*mm1m2m3 );
3329 c2dc4->Fill( pt, ms, tC2*tDC4, mm1*mpm3mqm1m2 );
3330 c4dc2->Fill( pt, ms, tC4*tDC2, mm1m2m3*mpmmq );
3331 c4dc4->Fill( pt, ms, tC4*tDC4, mm1m2m3*mpm3mqm1m2 );
3332 dc2dc4->Fill( pt, ms, tDC2*tDC4, mpmmq*mpm3mqm1m2 );
3335 // clean for next event
3336 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Reset();
3337 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Reset();
3338 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Reset();
3339 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Reset();
3340 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Reset();
3341 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Reset();
3343 //=======================================================================
3344 void AliAnalysisTaskFlowStrange::FillTrackVn(TString name,Double_t pt,Double_t phi,Double_t eta,Int_t fid) {
3346 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
3347 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
3348 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
3349 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
3351 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
3352 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
3353 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
3354 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
3355 Double_t qtpc = fQTPCA+fQTPCC;
3356 Double_t tpc_qmcos = (fQTPCACos+fQTPCCCos)/qtpc;
3357 Double_t tpc_qmsin = (fQTPCASin+fQTPCCSin)/qtpc;
3359 Double_t cosn = TMath::Cos(fHarmonic*phi);
3360 Double_t sinn = TMath::Sin(fHarmonic*phi);
3361 Double_t cos2n = TMath::Cos(2.0*fHarmonic*phi);
3362 Double_t sin2n = TMath::Sin(2.0*fHarmonic*phi);
3364 Double_t uQ, uQa, uQc, qaqc;
3365 // filling flow with vze
3366 qaqc = (vzea_qmcos*vzec_qmcos + vzea_qmsin*vzec_qmsin);
3367 uQa = (cosn*vzea_qmcos + sinn*vzea_qmsin);
3368 uQc = (cosn*vzec_qmcos + sinn*vzec_qmsin);
3369 Double_t cosmc = TMath::Cos( fHarmonic*GetMCDPHI(phi) );
3371 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQVZEA" ))->Fill( cosmc,uQa );
3372 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQVZEC" ))->Fill( cosmc,uQc );
3373 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI" ))->Fill( pt,cosmc );
3375 Double_t qaqt = (vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin);
3376 Double_t qcqt = (vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin);
3377 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" ))->Fill( pt,uQa,fQVZEA );
3378 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" ))->Fill( pt,uQc,fQVZEC );
3379 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" ))->Fill( pt,qaqc,fQVZEA*fQVZEC );
3380 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" ))->Fill( pt,qaqt,fQVZEA*qtpc );
3381 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" ))->Fill( pt,qcqt,fQVZEC*qtpc );
3383 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" ))->Fill( pt,uQa*uQc,fQVZEA*fQVZEC );
3384 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" ))->Fill( pt,uQa*qaqc,fQVZEA*fQVZEA*fQVZEC );
3385 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" ))->Fill( pt,uQc*qaqc,fQVZEC*fQVZEA*fQVZEC );
3386 // filling flow with tpc
3387 qaqc = (tpca_qmcos*tpcc_qmcos + tpca_qmsin*tpcc_qmsin);
3389 uQ = (cosn*tpca_qmcos + sinn*tpca_qmsin);
3390 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" ))->Fill( pt,uQ,fQTPCA );
3391 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" ))->Fill( pt,uQ*qaqc,fQTPCA*fQTPCA*fQTPCC );
3393 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQTPCA" ))->Fill( cosmc,uQ );
3396 uQ = (cosn*tpcc_qmcos + sinn*tpcc_qmsin);
3397 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" ))->Fill( pt,uQ,fQTPCC );
3398 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" ))->Fill( pt,uQ*qaqc,fQTPCC*fQTPCA*fQTPCC );
3400 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQTPCC" ))->Fill( cosmc,uQ );
3403 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" ))->Fill( pt,qaqc,fQTPCA*fQTPCC );
3405 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_P" ))->Fill( pt );
3406 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Fill( pt, cosn, 1.0 );
3407 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Fill( pt, sinn, 1.0 );
3409 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_Q" ))->Fill( pt );
3410 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Fill( pt, cosn, 1.0 );
3411 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Fill( pt, sinn, 1.0 );
3412 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Fill( pt, cos2n, 1.0 );
3413 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Fill( pt, sin2n, 1.0 );
3416 //=======================================================================
3417 void AliAnalysisTaskFlowStrange::FillDecayVn(TString name,Double_t ms,Double_t pt,Double_t phi,Double_t eta,Int_t fid1,Int_t fid2) {
3418 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "DecayYield_PtMass" ))->Fill( pt,ms );
3419 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "DecayAvgPt_PtMass" ))->Fill( pt,ms,pt );
3422 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
3423 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
3424 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
3425 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
3427 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
3428 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
3429 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
3430 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
3431 Double_t qtpc = fQTPCA+fQTPCC;
3432 Double_t tpc_qmcos = (fQTPCACos+fQTPCCCos)/qtpc;
3433 Double_t tpc_qmsin = (fQTPCASin+fQTPCCSin)/qtpc;
3435 Double_t cosn = TMath::Cos(fHarmonic*phi);
3436 Double_t sinn = TMath::Sin(fHarmonic*phi);
3437 Double_t cos2n = TMath::Cos(2.0*fHarmonic*phi);
3438 Double_t sin2n = TMath::Sin(2.0*fHarmonic*phi);
3440 Double_t uQ, uQa, uQc, qaqc;
3441 // filling flow with vze
3442 qaqc = (vzea_qmcos*vzec_qmcos + vzea_qmsin*vzec_qmsin);
3443 uQa = (cosn*vzea_qmcos + sinn*vzea_qmsin);
3444 uQc = (cosn*vzec_qmcos + sinn*vzec_qmsin);
3445 Double_t cosmc = TMath::Cos( fHarmonic*GetMCDPHI(phi) );
3447 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI" ))->Fill( pt,ms,cosmc );
3449 Double_t qaqt = (vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin);
3450 Double_t qcqt = (vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin);
3451 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" ))->Fill( pt,ms,uQa,fQVZEA );
3452 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" ))->Fill( pt,ms,uQc,fQVZEC );
3453 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" ))->Fill( pt,ms,qaqc,fQVZEA*fQVZEC );
3454 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" ))->Fill( pt,ms,qaqt,fQVZEA*qtpc );
3455 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" ))->Fill( pt,ms,qcqt,fQVZEC*qtpc );
3457 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" ))->Fill( pt,ms,uQa*uQc,fQVZEA*fQVZEC );
3458 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" ))->Fill( pt,ms,uQa*qaqc,fQVZEA*fQVZEA*fQVZEC );
3459 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" ))->Fill( pt,ms,uQc*qaqc,fQVZEC*fQVZEA*fQVZEC );
3460 // filling flow with tpc
3461 qaqc = (tpca_qmcos*tpcc_qmcos + tpca_qmsin*tpcc_qmsin);
3463 uQ = (cosn*tpca_qmcos + sinn*tpca_qmsin);
3464 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" ))->Fill( pt,ms,uQ,fQTPCA );
3465 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" ))->Fill( pt,ms,uQ*qaqc,fQTPCA*fQTPCA*fQTPCC );
3467 uQ = (cosn*tpcc_qmcos + sinn*tpcc_qmsin);
3468 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" ))->Fill( pt,ms,uQ,fQTPCC );
3469 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" ))->Fill( pt,ms,uQ*qaqc,fQTPCC*fQTPCA*fQTPCC );
3471 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" ))->Fill( pt,ms,qaqc,fQTPCA*fQTPCC );
3473 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_P" ))->Fill( pt,ms );
3474 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Fill( pt, ms, cosn, 1.0 );
3475 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Fill( pt, ms, sinn, 1.0 );
3476 if(InQTPC(fid1)||InQTPC(fid2)) {
3477 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_Q" ))->Fill( pt,ms );
3478 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Fill( pt, ms, cosn, 1.0 );
3479 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Fill( pt, ms, sinn, 1.0 );
3480 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Fill( pt, ms, cos2n, 1.0 );
3481 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Fill( pt, ms, sin2n, 1.0 );
3484 //=======================================================================
3485 Bool_t AliAnalysisTaskFlowStrange::InQTPC(Int_t id) {
3486 Bool_t ret = kFALSE;
3487 for(int i=0; i!=fQTPCA_nTracks; ++i)
3488 if(fQTPCA_fID[i]==id) {
3492 if(ret) return kTRUE;
3493 for(int i=0; i!=fQTPCC_nTracks; ++i)
3494 if(fQTPCC_fID[i]==id) {
3500 //=======================================================================
3501 void AliAnalysisTaskFlowStrange::ComputeTrackVn(TString name) {
3502 TProfile *uQa, *uQc, *qaqc, *uQaqaqc, *uQcqaqc;
3503 TArrayD *pasww, *pbsww, *pcsww;
3505 printf("<<%s>> SP TPC\n",name.Data());
3506 uQa = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" );
3507 uQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" );
3508 qaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" );
3509 uQaqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" );
3510 uQcqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" );
3511 pasww = uQa->GetBinSumw2();
3512 pbsww = uQc->GetBinSumw2();
3513 pcsww = qaqc->GetBinSumw2();
3515 TH1D *sptpca = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCA" );
3516 TH1D *sptpcc = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCC" );
3517 TH1D *sptpcaa = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCAA" );
3519 for(Int_t i=1; i!=sptpcaa->GetNbinsX()+1; ++i) {
3520 sptpcaa->SetBinContent(i,0);
3521 sptpcaa->SetBinError(i,0);
3522 sptpca->SetBinContent(i,0);
3523 sptpca->SetBinError(i,0);
3524 sptpcc->SetBinContent(i,0);
3525 sptpcc->SetBinError(i,0);
3526 double a = uQa->GetBinContent(i);
3527 double b = uQc->GetBinContent(i);
3528 double c = qaqc->GetBinContent(i);
3529 //if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3530 //if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3531 if(c<1e-100) continue;
3533 double vna = a/TMath::Sqrt(c);
3534 sptpca->SetBinContent(i,vna);
3536 double vnc = b/TMath::Sqrt(c);
3537 sptpcc->SetBinContent(i,vnc);
3539 double vn = (vna + vnc)/2.0;
3540 sptpcaa->SetBinContent(i,vn);
3542 double asw = uQa->GetBinEntries(i);
3543 double bsw = uQc->GetBinEntries(i);
3544 double csw = qaqc->GetBinEntries(i);
3545 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3546 double asww = pasww->At(i);
3547 double bsww = pbsww->At(i);
3548 double csww = pcsww->At(i);
3549 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3550 if((1<1e-100+asww/asw/asw)||(1<1e-100+bsww/bsw/bsw)||(1<1e-100+csww/csw/csw)) continue;
3551 if(TMath::AreEqualAbs(asww,asw*asw,1e-100)||
3552 TMath::AreEqualAbs(bsww,bsw*bsw,1e-100)||
3553 TMath::AreEqualAbs(csww,csw*csw,1e-100)) continue;
3554 double ac = uQaqaqc->GetBinContent(i);
3555 double bc = uQcqaqc->GetBinContent(i);
3556 double acsw = uQaqaqc->GetBinEntries(i);
3557 double bcsw = uQcqaqc->GetBinEntries(i);
3558 double ea = uQa->GetBinError(i)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3559 double eb = uQc->GetBinError(i)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3560 double ec = qaqc->GetBinError(i)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3561 //printf("%d >> ea^2 %.16f |||| asww %.16f | asw %.16f | 1-asww/asw/asw %.16f \n", i,ea, asww, asw, 1-asww/asw/asw);
3562 //printf("%d >> eb^2 %.16f |||| bsww %.16f | bsw %.16f | 1-bsww/bsw/bsw %.16f \n", i,eb, bsww, bsw, 1-bsww/bsw/bsw);
3563 //printf("%d >> ec^2 %.16f |||| csww %.16f | csw %.16f | 1-csww/csw/csw %.16f \n", i,ec, csww, csw, 1-csww/csw/csw);
3564 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3565 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3566 double evna = 1.0/TMath::Abs(c) * ( ea*ea + vna*vna/TMath::Abs(c)/4.0*ec*ec - vna/TMath::Sqrt(c)*eac );
3567 double evnc = 1.0/TMath::Abs(c) * ( eb*eb + vnc*vnc/TMath::Abs(c)/4.0*ec*ec - vnc/TMath::Sqrt(c)*ebc );
3568 //printf("%d >> evna^2 %.16f |||| ea %.16f | ec %.16f | eac %.16f | c %.16f\n", i,evna, ea, ec, eac,c);
3569 //printf("%d >> evnc^2 %.16f |||| eb %.16f | ec %.16f | ebc %.16f | c %.16f\n", i,evnc, eb, ec, ebc,c);
3570 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3571 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3572 sptpca->SetBinError(i,evna);
3573 sptpcc->SetBinError(i,evnc);
3574 sptpcaa->SetBinError(i,TMath::Sqrt(evna*evna+evnc*evnc)/2.0);
3577 printf("<<%s>> SP VZE\n",name.Data());
3578 double cvzea2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->GetBinContent( 1 );
3579 double cvzec2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->GetBinContent( 1 );
3580 if( TMath::AreEqualAbs(cvzea2+cvzec2,0,1e-100) ) return;
3581 uQa = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" );
3582 uQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" );
3583 qaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" );
3584 uQaqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" );
3585 uQcqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" );
3586 pasww = uQa->GetBinSumw2();
3587 pbsww = uQc->GetBinSumw2();
3588 pcsww = qaqc->GetBinSumw2();
3590 TProfile *qaqt = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" );
3591 TProfile *qcqt = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" );
3592 TProfile *uQauQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" );
3594 TH1D *spvzea = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEA" );
3595 TH1D *spvzec = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEC" );
3596 TH1D *spvzega = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEGA" );
3597 TH1D *spvzewa = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEWA" );
3598 for(Int_t i=1; i!=spvzewa->GetNbinsX()+1; ++i) {
3599 spvzega->SetBinContent(i,0);
3600 spvzega->SetBinError(i,0);
3601 spvzewa->SetBinContent(i,0);
3602 spvzewa->SetBinError(i,0);
3603 spvzea->SetBinContent(i,0);
3604 spvzea->SetBinError(i,0);
3605 spvzec->SetBinContent(i,0);
3606 spvzec->SetBinError(i,0);
3607 double asw = uQa->GetBinEntries(i);
3608 double bsw = uQc->GetBinEntries(i);
3609 double csw = qaqc->GetBinEntries(i);
3610 if(asw<1e-1||bsw<1e-1||csw<1e-1) continue;
3611 double asww = pasww->At(i);
3612 double bsww = pbsww->At(i);
3613 double csww = pcsww->At(i);
3614 if(asww<1e-1||bsww<1e-1||csww<1e-1) continue;
3615 if((1<asww/asw/asw)||(1<bsww/bsw/bsw)||(1<csww/csw/csw)) continue;
3616 double a = uQa->GetBinContent(i);
3617 double b = uQc->GetBinContent(i);
3618 double c = qaqc->GetBinContent(i);
3619 double at = qaqt->GetBinContent(i);
3620 double bt = qcqt->GetBinContent(i);
3621 if(TMath::AreEqualAbs(a,0,1e-10)) continue;
3622 if(TMath::AreEqualAbs(b,0,1e-10)) continue;
3623 if(TMath::AreEqualAbs(c,0,1e-10)) continue;
3624 if(TMath::AreEqualAbs(at,0,1e-10)) continue;
3625 if(TMath::AreEqualAbs(bt,0,1e-10)) continue;
3627 double aa = c*at/bt;
3628 if(aa<1e-100) continue;
3629 double vna = a/TMath::Sqrt(aa);
3630 spvzea->SetBinContent(i,vna);
3632 double bb = c*bt/at;
3633 if(bb<1e-100) continue;
3634 double vnc = b/TMath::Sqrt(bb);
3635 spvzec->SetBinContent(i,vnc);
3637 double vnwa = (cvzea2*vna + cvzec2*vnc) / (cvzea2+cvzec2);
3638 spvzewa->SetBinContent(i,vnwa);
3640 double vnga = a*b/c;
3641 if(vnga<1e-100) continue;
3642 vnga = TMath::Sqrt(vnga);
3643 spvzega->SetBinContent(i,vnga);
3645 double ab = uQauQc->GetBinContent(i);
3646 double ac = uQaqaqc->GetBinContent(i);
3647 double bc = uQcqaqc->GetBinContent(i);
3648 double absw = uQauQc->GetBinEntries(i);
3649 double acsw = uQaqaqc->GetBinEntries(i);
3650 double bcsw = uQcqaqc->GetBinEntries(i);
3651 double ea = uQa->GetBinError(i)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3652 double eb = uQc->GetBinError(i)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3653 double ec = qaqc->GetBinError(i)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3654 if(TMath::AreEqualAbs(1,absw/asw/bsw,1e-100)||TMath::AreEqualAbs(1,bcsw/bsw/csw,1e-100)||TMath::AreEqualAbs(1,acsw/asw/csw,1e-100)) continue;
3655 double eab = (ab-a*b)/(1-absw/asw/bsw)*absw/asw/bsw;
3656 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3657 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3658 double nc, nec, neac, nebc;
3663 double evna = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*ea*ea + a*a/4.0*nec*nec - a*TMath::Abs(nc)*neac*neac );
3668 double evnc = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*eb*eb + b*b/4.0*nec*nec - b*TMath::Abs(nc)*nebc*nebc );
3669 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3670 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3671 spvzea->SetBinError(i,evna);
3672 spvzec->SetBinError(i,evnc);
3673 double evnwa = TMath::Sqrt( cvzea2*cvzea2*evna*evna + cvzec2*cvzec2*evnc*evnc )/(cvzea2+cvzec2);
3674 spvzewa->SetBinError(i,evnwa);
3675 double evnga = 0.25/c/c * ( TMath::Abs(b*c/a)*ea*ea + TMath::Abs(a*c/b)*eb*eb + TMath::Abs(a*b/c)*ec*ec + 2*c*eab - 2*a*ebc - 2*b*eac );
3676 if(evnga>1e-100) evnga = TMath::Sqrt(evnga); else evnga=0;
3677 spvzega->SetBinError(i,evnga);
3679 printf("<<%s>> QC TPC\n",name.Data());
3681 TH1D *resC2 = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum2" );
3682 TH1D *resC4 = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum4" );
3683 TH1D *resDC2= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum2" );
3684 TH1D *resDC4= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum4" );
3685 TH1D *resvn2= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn2" );
3686 TH1D *resvn4= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn4" );
3688 TProfile *c2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
3689 TProfile *c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
3690 TProfile *dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
3691 TProfile *dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
3692 TProfile *c2c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
3693 TProfile *c2dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
3694 TProfile *c2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
3695 TProfile *c4dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
3696 TProfile *c4dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
3697 TProfile *dc2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
3698 TArrayD *c2sww = c2->GetBinSumw2();
3699 TArrayD *c4sww = c4->GetBinSumw2();
3700 TArrayD *dc2sww= dc2->GetBinSumw2();
3701 TArrayD *dc4sww= dc4->GetBinSumw2();
3702 for(Int_t i=1; i!=resvn2->GetNbinsX()+1; ++i) {
3704 double v_c2sw = c2->GetBinEntries(i);
3705 if(v_c2sw<1e-100) continue;
3706 double v_c2sww = c2sww->At(i);
3707 double v_c2 = c2->GetBinContent(i);
3708 double e_c2 = TMath::Sqrt(v_c2sww)/v_c2sw*c2->GetBinError(i);
3711 resC2->SetBinContent(i, cum2 );
3712 resC2->SetBinError(i, ecum2 );
3714 double v_c4sw = c4->GetBinEntries(i);
3715 if(v_c4sw<1e-100) continue;
3716 double v_c4sww = c4sww->At(i);
3717 double v_c4 = c4->GetBinContent(i);
3718 double e_c4 = TMath::Sqrt(v_c4sww)/v_c4sw*c4->GetBinError(i);
3719 double v_c2c4 = c2c4->GetBinContent(i);
3720 double v_c2c4sw = c2c4->GetBinEntries(i);
3721 if(TMath::AreEqualAbs(v_c2c4sw/v_c2sw/v_c4sw,1,1e-100)) continue;
3722 double covc2c4 = v_c2c4sw/v_c2sw/v_c4sw*(v_c2c4 - v_c2*v_c4)/(1-v_c2c4sw/v_c2sw/v_c4sw);
3723 double cum4 = v_c4 - 2*v_c2*v_c2;
3724 double ecum4= 16.0*v_c2*v_c2*e_c2*e_c2 + e_c4*e_c4 - 8.0*v_c2*covc2c4;
3725 if(ecum4<1e-100) continue;
3726 ecum4 = TMath::Sqrt( ecum4 );
3727 resC4->SetBinContent(i, cum4 );
3728 resC4->SetBinError(i, ecum4 );
3730 double v_dc2sw = dc2->GetBinEntries(i);
3731 if(v_dc2sw<1) continue;
3732 double v_dc2 = dc2->GetBinContent(i);
3733 double v_dc2sww = dc2sww->At(i);
3734 double e_dc2 = TMath::Sqrt(v_dc2sww)/v_dc2sw*dc2->GetBinError(i);
3735 double dcum2 = v_dc2;
3736 double edcum2= e_dc2;
3737 resDC2->SetBinContent(i, dcum2 );
3738 resDC2->SetBinError(i, edcum2 );
3740 if(v_c2<1e-100) continue;
3741 double dv22 = v_dc2/TMath::Sqrt(v_c2);
3742 double v_c2dc2 = c2dc2->GetBinContent(i);
3743 double v_c2dc2sw = c2dc2->GetBinEntries(i);
3744 if(TMath::AreEqualAbs(v_c2dc2sw/v_c2sw/v_dc2sw,1,1e-100)) continue;
3745 double covc2dc2 = v_c2dc2sw/v_c2sw/v_dc2sw*(v_c2dc2 - v_c2*v_dc2)/(1-v_c2dc2sw/v_c2sw/v_dc2sw);
3746 double edv22 = 0.25/v_c2/v_c2/v_c2*(v_dc2*v_dc2*e_c2*e_c2 + 4*v_c2*v_c2*e_dc2*e_dc2 - 4*v_c2*v_dc2*covc2dc2);
3747 //printf("%d >> dv22 %.16f || edv22^2 %.16f |||| v_c2dc2 %.16f | v_c2dc2sw %.16f | covc2dc2 %.16f \n", i,dv22,edv22,v_c2dc2,v_c2dc2sw,covc2dc2);
3748 if(edv22<1e-100) continue;
3749 edv22 = TMath::Sqrt(edv22);
3750 resvn2->SetBinContent(i,dv22);
3751 resvn2->SetBinError(i,edv22);
3753 double v_dc4sw = dc4->GetBinEntries(i);
3754 if(v_dc4sw<1) continue;
3755 double v_dc4 = dc4->GetBinContent(i);
3756 double v_dc4sww = dc4sww->At(i);
3757 double e_dc4 = TMath::Sqrt(v_dc4sww)/v_dc4sw*dc4->GetBinError(i);
3758 double dcum4 = v_dc4 - 2*v_c2*v_dc2;
3759 double v_c2dc4 = c2dc4->GetBinContent(i);
3760 double v_c2dc4sw = c2dc4->GetBinEntries(i);
3761 if(TMath::AreEqualAbs(v_c2dc4sw/v_c2sw/v_dc4sw,1,1e-100)) continue;
3762 double covc2dc4 = v_c2dc4sw/v_c2sw/v_dc4sw*(v_c2dc4 - v_c2*v_dc4)/(1-v_c2dc4sw/v_c2sw/v_dc4sw);
3763 double v_dc2dc4 = dc2dc4->GetBinContent(i);
3764 double v_dc2dc4sw = dc2dc4->GetBinEntries(i);
3765 if(TMath::AreEqualAbs(v_dc2dc4sw/v_dc2sw/v_dc4sw,1,1e-100)) continue;
3766 double covdc2dc4 = v_dc2dc4sw/v_dc2sw/v_dc4sw*(v_dc2dc4 - v_dc2*v_dc4)/(1-v_dc2dc4sw/v_dc2sw/v_dc4sw);
3767 double edcum4= ( +4.0*v_dc2*v_dc2*e_c2*e_c2
3768 +4.0*v_c2*v_c2*e_dc2*e_dc2
3770 +8.0*v_c2*v_dc2*covc2dc2
3772 -4.0*v_c2*covdc2dc4 );
3773 if(edcum4<1e-100) continue;
3774 edcum4 = TMath::Sqrt(edcum4);
3775 resDC4->SetBinContent(i, dcum4 );
3776 resDC4->SetBinError(i, edcum4 );
3778 if(cum4>1e-100) continue;
3779 double dv24 = -dcum4/TMath::Power(-cum4,0.75);
3780 double dterm1 = 2*v_c2*v_c2*v_dc2 - 3*v_c2*v_dc4 + 2*v_c4*v_dc2;
3781 double dterm2 = 9.0/16.0*dcum4*dcum4;
3782 double dterm3 = 4.0*v_c2*v_c2*cum4*cum4;
3783 double dterm4 = cum4*cum4;
3784 double dterm5 = -3.0/2.0*dcum4*dterm1;
3785 double dterm6 = -4.0*v_c2*cum4*dterm1;
3786 double dterm7 = -2.0*cum4*dterm1;
3787 double dterm8 = 3.0*v_c2*cum4*dcum4;
3788 double dterm9 = 3.0/2.0*cum4*dcum4;
3789 double dterm10= 4*v_c2*cum4*cum4;
3790 double v_c4dc2 = c4dc2->GetBinContent(i);
3791 double v_c4dc2sw = c4dc2->GetBinEntries(i);
3792 if(TMath::AreEqualAbs(v_c4dc2sw/v_c4sw/v_dc2sw,1,1e-100)) continue;
3793 double covc4dc2 = v_c4dc2sw/v_c4sw/v_dc2sw*(v_c4dc2 - v_c4*v_dc2)/(1-v_c4dc2sw/v_c4sw/v_dc2sw);
3794 double v_c4dc4 = c4dc4->GetBinContent(i);
3795 double v_c4dc4sw = c4dc4->GetBinEntries(i);
3796 if(TMath::AreEqualAbs(v_c4dc4sw/v_c4sw/v_dc4sw,1,1e-100)) continue;
3797 double covc4dc4 = v_c4dc4sw/v_c4sw/v_dc4sw*(v_c4dc4 - v_c4*v_dc4)/(1-v_c4dc4sw/v_c4sw/v_dc4sw);
3798 double edv24= 1.0/TMath::Power(-cum4,3.5)*(+dterm1*dterm1*e_c2*e_c2
3807 -dterm10*covdc2dc4);
3808 if(edv24<1e-100) continue;
3809 edv24 = TMath::Sqrt(edv24);
3810 resvn4->SetBinContent(i,dv24);
3811 resvn4->SetBinError(i,edv24);
3814 //=======================================================================
3815 void AliAnalysisTaskFlowStrange::ComputeDecayVn(TString name) {
3816 TProfile2D *uQa, *uQc, *qaqc, *uQaqaqc, *uQcqaqc;
3817 TArrayD *pasww, *pbsww, *pcsww;
3819 printf("<<%s>> SP TPC\n",name.Data());
3820 uQa = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" );
3821 uQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" );
3822 qaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" );
3823 uQaqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" );
3824 uQcqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" );
3825 pasww = uQa->GetBinSumw2();
3826 pbsww = uQc->GetBinSumw2();
3827 pcsww = qaqc->GetBinSumw2();
3829 TH2D *sptpca = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCA" );
3830 TH2D *sptpcc = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCC" );
3831 TH2D *sptpcaa = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCAA" );
3833 for(Int_t i=1; i!=sptpcaa->GetNbinsX()+1; ++i) {
3834 for(Int_t j=1; j!=sptpcaa->GetNbinsY()+1; ++j) {
3835 sptpcaa->SetBinContent(i,j,0);
3836 sptpcaa->SetBinError(i,j,0);
3837 sptpca->SetBinContent(i,j,0);
3838 sptpca->SetBinError(i,j,0);
3839 sptpcc->SetBinContent(i,j,0);
3840 sptpcc->SetBinError(i,j,0);
3841 double a = uQa->GetBinContent(i,j);
3842 double b = uQc->GetBinContent(i,j);
3843 double c = qaqc->GetBinContent(i,j);
3844 //if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3845 //if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3846 if(c<1e-100) {printf("skipping i=%d, j=%d due to c=%.16f\n",i,j,c); continue;}
3848 double vna = a/TMath::Sqrt(c);
3849 sptpca->SetBinContent(i,j,vna);
3851 double vnc = b/TMath::Sqrt(c);
3852 sptpcc->SetBinContent(i,j,vnc);
3854 double vn = (vna + vnc)/2.0;
3855 sptpcaa->SetBinContent(i,j,vn);
3857 int k = sptpcaa->GetBin(i,j);
3858 double asw = uQa->GetBinEntries(k);
3859 double bsw = uQc->GetBinEntries(k);
3860 double csw = qaqc->GetBinEntries(k);
3861 if(asw<1e-100||bsw<1e-100||csw<1e-100) {printf("skipping i=%d, j=%d due to asw=%f or bsw=%f or csw=%f\n",i,j,asw,bsw,csw); continue;}
3862 double asww = pasww->At(k);
3863 double bsww = pbsww->At(k);
3864 double csww = pcsww->At(k);
3865 if(asww<1e-100||bsww<1e-100||csww<1e-100) {printf("skipping i=%d, j=%d due to asww=%f or bsww=%f or csww=%f\n",i,j,asww,bsww,csww); continue;}
3866 if((1<1e-100+asww/asw/asw)||(1<1e-100+bsww/bsw/bsw)||(1<1e-100+csww/csw/csw)) {printf("skipping i=%d, j=%d due to COVa=%f or COVb=%f or COVc=%f\n",i,j,asww/asw/asw,bsww/bsw/bsw,csww/csw/csw); continue;}
3867 if(TMath::AreEqualAbs(asww,asw*asw,1e-100)||
3868 TMath::AreEqualAbs(bsww,bsw*bsw,1e-100)||
3869 TMath::AreEqualAbs(csww,csw*csw,1e-100)) {printf("skipping i=%d, j=%d due to funny coincidence\n",i,j); continue;}
3870 double ac = uQaqaqc->GetBinContent(i,j);
3871 double bc = uQcqaqc->GetBinContent(i,j);
3872 double acsw = uQaqaqc->GetBinEntries(k);
3873 double bcsw = uQcqaqc->GetBinEntries(k);
3874 double ea = uQa->GetBinError(i,j)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3875 double eb = uQc->GetBinError(i,j)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3876 double ec = qaqc->GetBinError(i,j)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3877 //printf("%d >> ea^2 %.16f |||| asww %.16f | asw %.16f | 1-asww/asw/asw %.16f \n", i,ea, asww, asw, 1-asww/asw/asw);
3878 //printf("%d >> eb^2 %.16f |||| bsww %.16f | bsw %.16f | 1-bsww/bsw/bsw %.16f \n", i,eb, bsww, bsw, 1-bsww/bsw/bsw);
3879 //printf("%d >> ec^2 %.16f |||| csww %.16f | csw %.16f | 1-csww/csw/csw %.16f \n", i,ec, csww, csw, 1-csww/csw/csw);
3880 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3881 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3882 double evna = 1.0/TMath::Abs(c) * ( ea*ea + vna*vna/TMath::Abs(c)/4.0*ec*ec - vna/TMath::Sqrt(c)*eac );
3883 double evnc = 1.0/TMath::Abs(c) * ( eb*eb + vnc*vnc/TMath::Abs(c)/4.0*ec*ec - vnc/TMath::Sqrt(c)*ebc );
3884 //printf("%d >> evna^2 %.16f |||| ea %.16f | ec %.16f | eac %.16f | c %.16f\n", i,evna, ea, ec, eac,c);
3885 //printf("%d >> evnc^2 %.16f |||| eb %.16f | ec %.16f | ebc %.16f | c %.16f\n", i,evnc, eb, ec, ebc,c);
3886 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3887 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3888 sptpca->SetBinError(i,j,evna);
3889 sptpcc->SetBinError(i,j,evnc);
3890 sptpcaa->SetBinError(i,j,TMath::Sqrt(evna*evna+evnc*evnc)/2.0);
3894 printf("<<%s>> SP VZE\n",name.Data());
3895 double cvzea2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->GetBinContent( 1 );
3896 double cvzec2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->GetBinContent( 1 );
3897 if( TMath::AreEqualAbs(cvzea2+cvzec2,0,1e-100) ) return;
3898 uQa = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" );
3899 uQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" );
3900 qaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" );
3901 uQaqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" );
3902 uQcqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" );
3903 pasww = uQa->GetBinSumw2();
3904 pbsww = uQc->GetBinSumw2();
3905 pcsww = qaqc->GetBinSumw2();
3907 TProfile2D *qaqt = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" );
3908 TProfile2D *qcqt = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" );
3909 TProfile2D *uQauQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" );
3911 TH2D *spvzea = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEA" );
3912 TH2D *spvzec = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEC" );
3913 TH2D *spvzega = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEGA" );
3914 TH2D *spvzewa = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEWA" );
3915 for(Int_t i=1; i!=spvzewa->GetNbinsX()+1; ++i) {
3916 for(Int_t j=1; j!=spvzewa->GetNbinsY()+1; ++j) {
3917 spvzega->SetBinContent(i,j,0);
3918 spvzega->SetBinError(i,j,0);
3919 spvzewa->SetBinContent(i,j,0);
3920 spvzewa->SetBinError(i,j,0);
3921 spvzea->SetBinContent(i,j,0);
3922 spvzea->SetBinError(i,j,0);
3923 spvzec->SetBinContent(i,j,0);
3924 spvzec->SetBinError(i,j,0);
3925 double a = uQa->GetBinContent(i,j);
3926 double b = uQc->GetBinContent(i,j);
3927 double c = qaqc->GetBinContent(i,j);
3928 double at = qaqt->GetBinContent(i,j);
3929 double bt = qcqt->GetBinContent(i,j);
3930 if(TMath::AreEqualAbs(a,0,1e-100)) {printf("skipping A\n"); continue;}
3931 if(TMath::AreEqualAbs(b,0,1e-100)) {printf("skipping B\n"); continue;}
3932 if(TMath::AreEqualAbs(c,0,1e-100)) {printf("skipping C\n"); continue;}
3933 if(TMath::AreEqualAbs(at,0,1e-100)) {printf("skipping AT\n"); continue;}
3934 if(TMath::AreEqualAbs(bt,0,1e-100)) {printf("skipping CT\n"); continue;}
3936 double aa = c*at/bt;
3937 if(aa<1e-100) {printf("AA\n"); continue;}
3938 double vna = a/TMath::Sqrt(aa);
3939 spvzea->SetBinContent(i,j,vna);
3941 double bb = c*bt/at;
3942 if(bb<1e-100) {printf("BB\n"); continue;}
3943 double vnc = b/TMath::Sqrt(bb);
3944 spvzec->SetBinContent(i,j,vnc);
3946 double vnwa = (cvzea2*vna + cvzec2*vnc) / (cvzea2+cvzec2);
3947 spvzewa->SetBinContent(i,j,vnwa);
3949 double vnga = a*b/c;
3950 if(vnga<1e-100) continue;
3951 vnga = TMath::Sqrt(vnga);
3952 spvzega->SetBinContent(i,j,vnga);
3954 int k = spvzea->GetBin(i,j);
3955 double asw = uQa->GetBinEntries(k);
3956 double bsw = uQc->GetBinEntries(k);
3957 double csw = qaqc->GetBinEntries(k);
3958 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3959 double asww = pasww->At(k);
3960 double bsww = pbsww->At(k);
3961 double csww = pcsww->At(k);
3962 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3963 if((1<asww/asw/asw)||(1<bsww/bsw/bsw)||(1<csww/csw/csw)) continue;
3964 double ab = uQauQc->GetBinContent(i,j);
3965 double ac = uQaqaqc->GetBinContent(i,j);
3966 double bc = uQcqaqc->GetBinContent(i,j);
3967 double absw = uQauQc->GetBinEntries(k);
3968 double acsw = uQaqaqc->GetBinEntries(k);
3969 double bcsw = uQcqaqc->GetBinEntries(k);
3970 if(TMath::AreEqualAbs(1,absw/asw/bsw,1e-100)||TMath::AreEqualAbs(1,bcsw/bsw/csw,1e-100)||TMath::AreEqualAbs(1,acsw/asw/csw,1e-100)) continue;
3971 double ea = uQa->GetBinError(i,j)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3972 double eb = uQc->GetBinError(i,j)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3973 double ec = qaqc->GetBinError(i,j)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3974 double eab = (ab-a*b)/(1-absw/asw/bsw)*absw/asw/bsw;
3975 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3976 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3977 double nc, nec, neac, nebc;
3982 double evna = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*ea*ea + a*a/4.0*nec*nec - a*TMath::Abs(nc)*neac*neac );
3987 double evnc = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*eb*eb + b*b/4.0*nec*nec - b*TMath::Abs(nc)*nebc*nebc );
3988 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3989 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3990 spvzea->SetBinError(i,j,evna);
3991 spvzec->SetBinError(i,j,evnc);
3992 double evnwa = TMath::Sqrt( cvzea2*cvzea2*evna*evna + cvzec2*cvzec2*evnc*evnc )/(cvzea2+cvzec2);
3993 spvzewa->SetBinError(i,j,evnwa);
3994 double evnga = 0.25/c/c * ( TMath::Abs(b*c/a)*ea*ea + TMath::Abs(a*c/b)*eb*eb + TMath::Abs(a*b/c)*ec*ec + 2*c*eab - 2*a*ebc - 2*b*eac );
3995 if(evnga>1e-100) evnga = TMath::Sqrt(evnga); else evnga=0;
3996 spvzega->SetBinError(i,j,evnga);
3999 printf("<<%s>> QC TPC\n",name.Data());
4001 TH2D *resC2 = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum2" );
4002 TH2D *resC4 = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum4" );
4003 TH2D *resDC2= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum2" );
4004 TH2D *resDC4= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum4" );
4005 TH2D *resvn2= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn2" );
4006 TH2D *resvn4= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn4" );
4008 TProfile2D *c2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
4009 TProfile2D *c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
4010 TProfile2D *dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
4011 TProfile2D *dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
4012 TProfile2D *c2c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
4013 TProfile2D *c2dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
4014 TProfile2D *c2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
4015 TProfile2D *c4dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
4016 TProfile2D *c4dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
4017 TProfile2D *dc2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
4018 TArrayD *c2sww = c2->GetBinSumw2();
4019 TArrayD *c4sww = c4->GetBinSumw2();
4020 TArrayD *dc2sww= dc2->GetBinSumw2();
4021 TArrayD *dc4sww= dc4->GetBinSumw2();
4022 for(Int_t i=1; i!=resvn2->GetNbinsX()+1; ++i) {
4023 for(Int_t j=1; j!=resvn2->GetNbinsY()+1; ++j) {
4025 int k = c2->GetBin(i,j);
4026 double v_c2sw = c2->GetBinEntries(k);
4027 if(v_c2sw<1e-100) continue;
4028 double v_c2sww = c2sww->At(k);
4029 double v_c2 = c2->GetBinContent(i,j);
4030 double e_c2 = TMath::Sqrt(v_c2sww)/v_c2sw*c2->GetBinError(i,j);
4033 resC2->SetBinContent(i,j, cum2 );
4034 resC2->SetBinError(i,j, ecum2 );
4036 double v_c4sw = c4->GetBinEntries(k);
4037 if(v_c4sw<1e-100) continue;
4038 double v_c4sww = c4sww->At(k);
4039 double v_c4 = c4->GetBinContent(i,j);
4040 double e_c4 = TMath::Sqrt(v_c4sww)/v_c4sw*c4->GetBinError(i,j);
4041 double v_c2c4 = c2c4->GetBinContent(i,j);
4042 double v_c2c4sw = c2c4->GetBinEntries(k);
4043 if(TMath::AreEqualAbs(v_c2c4sw/v_c2sw/v_c4sw,1,1e-100)) continue;
4044 double covc2c4 = v_c2c4sw/v_c2sw/v_c4sw*(v_c2c4 - v_c2*v_c4)/(1-v_c2c4sw/v_c2sw/v_c4sw);
4045 double cum4 = v_c4 - 2*v_c2*v_c2;
4046 double ecum4= 16.0*v_c2*v_c2*e_c2*e_c2 + e_c4*e_c4 - 8.0*v_c2*covc2c4;
4047 if(ecum4<1e-100) continue;
4048 ecum4 = TMath::Sqrt( ecum4 );
4049 resC4->SetBinContent(i,j, cum4 );
4050 resC4->SetBinError(i,j, ecum4 );
4052 double v_dc2sw = dc2->GetBinEntries(k);
4053 if(v_dc2sw<1) continue;
4054 double v_dc2 = dc2->GetBinContent(i,j);
4055 double v_dc2sww = dc2sww->At(k);
4056 double e_dc2 = TMath::Sqrt(v_dc2sww)/v_dc2sw*dc2->GetBinError(i,j);
4057 double dcum2 = v_dc2;
4058 double edcum2= e_dc2;
4059 resDC2->SetBinContent(i,j, dcum2 );
4060 resDC2->SetBinError(i,j, edcum2 );
4062 if(v_c2<1e-100) continue;
4063 double dv22 = v_dc2/TMath::Sqrt(v_c2);
4064 double v_c2dc2 = c2dc2->GetBinContent(i,j);
4065 double v_c2dc2sw = c2dc2->GetBinEntries(k);
4066 if(TMath::AreEqualAbs(v_c2dc2sw/v_c2sw/v_dc2sw,1,1e-100)) continue;
4067 double covc2dc2 = v_c2dc2sw/v_c2sw/v_dc2sw*(v_c2dc2 - v_c2*v_dc2)/(1-v_c2dc2sw/v_c2sw/v_dc2sw);
4068 double edv22 = 0.25/v_c2/v_c2/v_c2*(v_dc2*v_dc2*e_c2*e_c2 + 4*v_c2*v_c2*e_dc2*e_dc2 - 4*v_c2*v_dc2*covc2dc2);
4069 //printf("%d >> dv22 %.16f || edv22^2 %.16f |||| v_c2dc2 %.16f | v_c2dc2sw %.16f | covc2dc2 %.16f \n", i,dv22,edv22,v_c2dc2,v_c2dc2sw,covc2dc2);
4070 if(edv22<1e-100) continue;
4071 edv22 = TMath::Sqrt(edv22);
4072 resvn2->SetBinContent(i,j,dv22);
4073 resvn2->SetBinError(i,j,edv22);
4075 double v_dc4sw = dc4->GetBinEntries(k);
4076 if(v_dc4sw<1) continue;
4077 double v_dc4 = dc4->GetBinContent(i,j);
4078 double v_dc4sww = dc4sww->At(k);
4079 double e_dc4 = TMath::Sqrt(v_dc4sww)/v_dc4sw*dc4->GetBinError(i,j);
4080 double dcum4 = v_dc4 - 2*v_c2*v_dc2;
4081 double v_c2dc4 = c2dc4->GetBinContent(i,j);
4082 double v_c2dc4sw = c2dc4->GetBinEntries(k);
4083 if(TMath::AreEqualAbs(v_c2dc4sw/v_c2sw/v_dc4sw,1,1e-100)) continue;
4084 double covc2dc4 = v_c2dc4sw/v_c2sw/v_dc4sw*(v_c2dc4 - v_c2*v_dc4)/(1-v_c2dc4sw/v_c2sw/v_dc4sw);
4085 double v_dc2dc4 = dc2dc4->GetBinContent(i,j);
4086 double v_dc2dc4sw = dc2dc4->GetBinEntries(k);
4087 if(TMath::AreEqualAbs(v_dc2dc4sw/v_dc2sw/v_dc4sw,1,1e-100)) continue;
4088 double covdc2dc4 = v_dc2dc4sw/v_dc2sw/v_dc4sw*(v_dc2dc4 - v_dc2*v_dc4)/(1-v_dc2dc4sw/v_dc2sw/v_dc4sw);
4089 double edcum4= ( +4.0*v_dc2*v_dc2*e_c2*e_c2
4090 +4.0*v_c2*v_c2*e_dc2*e_dc2
4092 +8.0*v_c2*v_dc2*covc2dc2
4094 -4.0*v_c2*covdc2dc4 );
4095 if(edcum4<1e-100) continue;
4096 edcum4 = TMath::Sqrt(edcum4);
4097 resDC4->SetBinContent(i,j, dcum4 );
4098 resDC4->SetBinError(i,j, edcum4 );
4100 if(cum4>1e-100) continue;
4101 double dv24 = -dcum4/TMath::Power(-cum4,0.75);
4102 double dterm1 = 2*v_c2*v_c2*v_dc2 - 3*v_c2*v_dc4 + 2*v_c4*v_dc2;
4103 double dterm2 = 9.0/16.0*dcum4*dcum4;
4104 double dterm3 = 4.0*v_c2*v_c2*cum4*cum4;
4105 double dterm4 = cum4*cum4;
4106 double dterm5 = -3.0/2.0*dcum4*dterm1;
4107 double dterm6 = -4.0*v_c2*cum4*dterm1;
4108 double dterm7 = -2.0*cum4*dterm1;
4109 double dterm8 = 3.0*v_c2*cum4*dcum4;
4110 double dterm9 = 3.0/2.0*cum4*dcum4;
4111 double dterm10= 4*v_c2*cum4*cum4;
4112 double v_c4dc2 = c4dc2->GetBinContent(i,j);
4113 double v_c4dc2sw = c4dc2->GetBinEntries(k);
4114 if(TMath::AreEqualAbs(v_c4dc2sw/v_c4sw/v_dc2sw,1,1e-100)) continue;
4115 double covc4dc2 = v_c4dc2sw/v_c4sw/v_dc2sw*(v_c4dc2 - v_c4*v_dc2)/(1-v_c4dc2sw/v_c4sw/v_dc2sw);
4116 double v_c4dc4 = c4dc4->GetBinContent(i,j);
4117 double v_c4dc4sw = c4dc4->GetBinEntries(k);
4118 if(TMath::AreEqualAbs(v_c4dc4sw/v_c4sw/v_dc4sw,1,1e-100)) continue;
4119 double covc4dc4 = v_c4dc4sw/v_c4sw/v_dc4sw*(v_c4dc4 - v_c4*v_dc4)/(1-v_c4dc4sw/v_c4sw/v_dc4sw);
4120 double edv24= 1.0/TMath::Power(-cum4,3.5)*(+dterm1*dterm1*e_c2*e_c2
4129 -dterm10*covdc2dc4);
4130 if(edv24<1e-100) continue;
4131 edv24 = TMath::Sqrt(edv24);
4132 resvn4->SetBinContent(i,j,dv24);
4133 resvn4->SetBinError(i,j,edv24);
4138 //=======================================================================
4139 void AliAnalysisTaskFlowStrange::OpenToyModel() {
4140 fList = new TList();
4144 tList=new TList(); tList->SetName("ToyVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
4152 //=======================================================================
4153 void AliAnalysisTaskFlowStrange::CloseToyModel() {
4155 ComputeDecayVn("ToyVn");
4157 //=======================================================================
4158 void AliAnalysisTaskFlowStrange::MakeToyEvent( Int_t seed, Int_t m_decay,Double_t v_decay,
4159 Double_t mass_decay_mu,Double_t mass_decay_sg,
4160 Int_t m_bgr,Double_t v_bgr,
4161 Int_t mtpc_a,Double_t v_tpca,Int_t mtpc_c,Double_t v_tpcc,
4162 Int_t mvze_a,Double_t v_vzea,Int_t mvze_c,Double_t v_vzec ) {
4163 gRandom->SetSeed( seed );
4165 fMCEP = gRandom->Rndm()*TMath::Pi();
4166 TF1 tf1_tpca( "dphitpca", Form("1+2*%f*TMath::Cos(2*x)",v_tpca),0,TMath::TwoPi() );
4167 TF1 tf1_tpcc( "dphitpcc", Form("1+2*%f*TMath::Cos(2*x)",v_tpcc),0,TMath::TwoPi() );
4168 TF1 tf1_vzea( "dphivzea", Form("1+2*%f*TMath::Cos(2*x)",v_vzea),0,TMath::TwoPi() );
4169 TF1 tf1_vzec( "dphivzec", Form("1+2*%f*TMath::Cos(2*x)",v_vzec),0,TMath::TwoPi() );
4170 TF1 tf1_decay( "dphidecay", Form("1+2*%f*TMath::Cos(2*x)",v_decay),0,TMath::TwoPi() );
4171 TF1 tf1_bgr( "dphibgr", Form("1+2*%f*TMath::Cos(2*x)",v_bgr),0,TMath::TwoPi() );
4173 fQTPCACos=fQTPCASin=fQTPCA=0;
4174 fQTPCCCos=fQTPCCSin=fQTPCC=0;
4175 fQTPC2hCos=fQTPC2hSin=0;
4176 fQVZEACos=fQVZEASin=fQVZEA=0;
4177 fQVZECCos=fQVZECSin=fQVZEC=0;
4178 for(int m=0; m!=mtpc_a; ++m) {
4179 phi = tf1_tpca.GetRandom() + fMCEP;
4180 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4181 eta = gRandom->Rndm()*(fRFPAmaxEta-fRFPAminEta)+fRFPAminEta;
4182 fQTPCACos += TMath::Cos(fHarmonic*phi);
4183 fQTPCASin += TMath::Sin(fHarmonic*phi);
4185 fQTPC2hCos += TMath::Cos(2*fHarmonic*phi);
4186 fQTPC2hSin += TMath::Sin(2*fHarmonic*phi);
4187 fQTPCA_fID[m] = -99;
4188 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta );
4190 for(int m=0; m!=mtpc_c; ++m) {
4191 phi = tf1_tpcc.GetRandom() + fMCEP;
4192 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4193 eta = gRandom->Rndm()*(fRFPCmaxEta-fRFPCminEta)+fRFPCminEta;
4194 fQTPCCCos += TMath::Cos(fHarmonic*phi);
4195 fQTPCCSin += TMath::Sin(fHarmonic*phi);
4197 fQTPC2hCos += TMath::Cos(2*fHarmonic*phi);
4198 fQTPC2hSin += TMath::Sin(2*fHarmonic*phi);
4199 fQTPCC_fID[m] = -99;
4200 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta );
4202 for(int m=0; m!=mvze_a; ++m) {
4203 phi = tf1_vzea.GetRandom() + fMCEP;
4204 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4205 eta = gRandom->Rndm()*2-3.5;
4206 fQVZEACos += TMath::Cos(fHarmonic*phi);
4207 fQVZEASin += TMath::Sin(fHarmonic*phi);
4209 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta );
4211 for(int m=0; m!=mvze_c; ++m) {
4212 phi = tf1_vzec.GetRandom() + fMCEP;
4213 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4214 eta = gRandom->Rndm()*2+2.5;
4215 fQVZECCos += TMath::Cos(fHarmonic*phi);
4216 fQVZECSin += TMath::Sin(fHarmonic*phi);
4218 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta );
4220 fQTPCA_nTracks = mtpc_a;
4221 fQTPCC_nTracks = mtpc_c;
4225 double ptrange = fPtBinEdge[fPtBins] - fPtBinEdge[0];
4227 for(int m=0; m!=m_decay; ++m) {
4228 phi = tf1_decay.GetRandom() + fMCEP;
4229 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4230 eta = gRandom->Rndm()*1.6-0.8;
4231 pt = gRandom->Rndm()*ptrange + fPtBinEdge[0];
4232 mass = gRandom->Gaus(mass_decay_mu,mass_decay_sg);
4233 FillDecayVn("ToyVn",mass,pt,phi,eta,+999,+999);
4235 for(int m=0; m!=m_bgr; ++m) {
4236 phi = tf1_bgr.GetRandom() + fMCEP;
4237 if(phi>TMath::TwoPi()) phi -= TMath::TwoPi();
4238 eta = gRandom->Rndm()*1.6-0.8;
4239 pt = gRandom->Rndm()*ptrange + fPtBinEdge[0];
4240 mass = gRandom->Rndm()*(fMaxMass-fMinMass)+fMinMass;
4241 FillDecayVn("ToyVn",mass,pt,phi,eta,+999,+999);
4243 QCStoreDecayVn("ToyVn");