-modified qa task
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / macrosJPSI / AddTask_jpsi_Default.C
1 void InitHistograms(AliDielectron *die, Int_t cutDefinition);
2 void InitCF(AliDielectron* die, Int_t cutDefinition);
3
4 void SetupEventCuts(AliDielectron *die, Int_t cutDefinition);
5 void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition);
6 void SetupPairCuts(AliDielectron *die,  Int_t cutDefinition);
7
8 void SetupMCsignals(AliDielectron *die);
9 TVectorD *GetRunNumbers();
10
11 TString names=("default");
12 TObjArray *arrNames=names.Tokenize(";");
13
14 const Int_t nDie=arrNames->GetEntries();
15 Bool_t isAOD=kFALSE;
16 Bool_t hasMC=kFALSE;
17 Int_t iPeriod=-1;
18 enum { k10b=0, k10c, k10d, k10e, k10f, k10h, k11a, k11d, k11h, k12h };
19
20 //______________________________________________________________________________________
21 AliAnalysisTask* AddTask_jpsi_Default(TString prod="", Bool_t isMC=kFALSE)
22 {
23   //get the current analysis manager
24
25   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
26   if (!mgr) {
27     Error("AddTask_jpsi_JPsi", "No analysis manager found.");
28     return 0;
29   }
30
31   //Do we have an MC handler?
32   hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0);
33
34   //AOD input?
35   isAOD=mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
36   if(isAOD) hasMC=isMC;
37
38   //Get the current train configuration
39   TString trainConfig=gSystem->Getenv("CONFIG_FILE");
40   TString list=gSystem->Getenv("LIST");
41   if( list.IsNull()) list=prod;
42
43   // selected period
44   if(      !prod.CompareTo("LHC10b") ) iPeriod = k10b;
45   else if( !prod.CompareTo("LHC10c") ) iPeriod = k10c;
46   else if( !prod.CompareTo("LHC10d") ) iPeriod = k10d;
47   else if( !prod.CompareTo("LHC10e") ) iPeriod = k10e;
48   else if( !prod.CompareTo("LHC10f") ) iPeriod = k10f;
49   else if( !prod.CompareTo("LHC10h") ) iPeriod = k10h;
50   else if( !prod.CompareTo("LHC11a") ) iPeriod = k11a;
51   else if( !prod.CompareTo("LHC11d") ) iPeriod = k11d;
52   else if( !prod.CompareTo("LHC11h") ) iPeriod = k11h;
53   else if( !prod.CompareTo("LHC12h") ) iPeriod = k12h;
54
55   // // aod monte carlo
56   // if( list.Contains("LHC11a10") ||
57   //     list.Contains("LHC11b10") ||
58   //     list.Contains("LHC12a17") ||
59   //     list.Contains("fix")
60   //     ) hasMC=kTRUE;
61
62   //create task and add it to the manager
63   AliAnalysisTaskMultiDielectron *task=new AliAnalysisTaskMultiDielectron("JpsiDefault");
64   //  task->SetBeamEnergy(1380.); // not neeeded since we are not looking at helicity and Collins-Soper coordinates
65   if (!hasMC) task->UsePhysicsSelection();
66
67   // add special triggers
68   switch(iPeriod) {
69   case k11d: task->SetTriggerMask(AliVEvent::kEMCEJE+AliVEvent::kEMC7+AliVEvent::kEMCEGA);     break;
70   case k11h: task->SetTriggerMask(AliVEvent::kMB+AliVEvent::kCentral+AliVEvent::kSemiCentral); break;
71   case k12h: task->SetTriggerMask(AliVEvent::kAnyINT);                                         break;
72   }
73   mgr->AddTask(task);
74
75   //add dielectron analysis with different cuts to the task
76   for (Int_t i=0; i<nDie; ++i){ //nDie defined in config file
77     AliDielectron *jpsi=ConfigDefault(i);
78     if (!jpsi) continue;
79     jpsi->SetHasMC(hasMC);
80     task->AddDielectron(jpsi);
81   }
82
83   //   task->SetTriggerOnV0AND();
84   //   if ( trainConfig=="pp" ) task->SetRejectPileup();
85   
86   //create output container
87   AliAnalysisDataContainer *coutput1 =
88     mgr->CreateContainer("jpsi_Default_tree",
89                          TTree::Class(),
90                          AliAnalysisManager::kExchangeContainer,
91                          "jpsi_Default_default");
92   
93   AliAnalysisDataContainer *cOutputHist1 =
94     mgr->CreateContainer("jpsi_Default_QA",
95                          TList::Class(),
96                          AliAnalysisManager::kOutputContainer,
97                          "jpsi_Default.root");
98
99   AliAnalysisDataContainer *cOutputHist2 =
100     mgr->CreateContainer("jpsi_Default_CF",
101                          TList::Class(),
102                          AliAnalysisManager::kOutputContainer,
103                          "jpsi_Default.root");
104   
105   AliAnalysisDataContainer *cOutputHist3 =
106     mgr->CreateContainer("jpsi_Default_EventStat",
107                          TH1D::Class(),
108                          AliAnalysisManager::kOutputContainer,
109                          "jpsi_Default.root");
110   
111   mgr->ConnectInput(task,  0, mgr->GetCommonInputContainer());
112   mgr->ConnectOutput(task, 0, coutput1 );
113   mgr->ConnectOutput(task, 1, cOutputHist1);
114   mgr->ConnectOutput(task, 2, cOutputHist2);
115   mgr->ConnectOutput(task, 3, cOutputHist3);
116   
117   return task;
118 }
119
120
121 //______________________________________________________________________________________
122 //______________________________________________________________________________________
123 //______________________________________________________________________________________
124 //
125 // Here the configuration part starts
126 //
127 AliDielectron* ConfigDefault(Int_t cutDefinition)
128 {
129   //
130   // Setup the instance of AliDielectron
131   //
132   
133   // create the actual framework object
134   TString name=Form("%02d",cutDefinition);
135   if (cutDefinition<arrNames->GetEntriesFast()){
136     name=arrNames->At(cutDefinition)->GetName();
137   }
138   AliDielectron *die =
139   new AliDielectron(Form("%s",name.Data()),
140                     Form("Track cuts: %s",name.Data()));
141   
142   /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
143   SetupEventCuts(die,cutDefinition);
144   SetupTrackCuts(die,cutDefinition);
145   SetupPairCuts(die,cutDefinition);
146
147   /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv MISC vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
148   // Monte Carlo Signals
149   if (hasMC) SetupMCsignals(die);
150   // prefilter settings
151   die->SetPreFilterUnlikeOnly();//  die->SetNoPairing();//  die->SetPreFilterAllSigns();
152   // cut QA
153   die->SetCutQA();
154
155   /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv OUTPUT vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
156   InitHistograms(die,cutDefinition);
157   InitCF(die,cutDefinition);
158
159
160   //   AliDielectronMixingHandler *mix=new AliDielectronMixingHandler;
161   //   mix->AddVariable(AliDielectronVarManager::kZvPrim,"-10,-8,-5,0,5,8,10");
162   //   mix->AddVariable(AliDielectronVarManager::kPhi ,"0,3.14,6.3");
163   //   mix->SetDepth(10);
164   //  die->SetMixingHandler(mix);
165   //
166   
167   return die;
168 }
169
170 //______________________________________________________________________________________
171 void SetupEventCuts(AliDielectron *die, Int_t cutDefinition)
172 {
173   //
174   // Setup the event cuts
175   //
176
177   AliDielectronEventCuts *eventCuts=new AliDielectronEventCuts("eventCuts","eventCuts");
178   if(isAOD) eventCuts->SetVertexType(AliDielectronEventCuts::kVtxAny);
179   eventCuts->SetRequireVertex();
180   eventCuts->SetMinVtxContributors(1);
181   eventCuts->SetVertexZ(-10.,+10.);
182   eventCuts->Print();
183   die->GetEventFilter().AddCuts(eventCuts);
184
185 }
186
187 //______________________________________________________________________________________
188 void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition)
189 {
190   //
191   // Setup the track cuts
192   //
193
194   /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv TRACK CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
195   AliDielectronVarCuts *varAccCuts   = new AliDielectronVarCuts("acc","acc");
196   varAccCuts->AddCut(AliDielectronVarManager::kPt,           0.8, 1e30);
197   varAccCuts->AddCut(AliDielectronVarManager::kEta,         -0.9,  0.9);
198   die->GetTrackFilter().AddCuts(varAccCuts);
199   varAccCuts->Print();
200
201
202   AliDielectronVarCuts   *varRecCuts = new AliDielectronVarCuts("VarRecCuts","VarRecCuts");
203   varRecCuts->AddCut(AliDielectronVarManager::kNclsTPC,      70.,   160.);
204   varRecCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
205   varRecCuts->AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
206   varRecCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
207   //  varRecCuts->AddCut(AliDielectronVarManager::kITSchi2Cl,     0.  ,   36.     ); // not defined in AOD
208   varRecCuts->AddCut(AliDielectronVarManager::kKinkIndex0,   0.0);
209
210   AliDielectronTrackCuts *trkRecCuts = new AliDielectronTrackCuts("TrkRecCuts","TrkRecCuts");
211   trkRecCuts->SetITSclusterCut(AliDielectronTrackCuts::kOneOf, 3); // SPD any
212   trkRecCuts->SetRequireITSRefit(kTRUE);
213   trkRecCuts->SetRequireTPCRefit(kTRUE);
214
215   AliDielectronCutGroup  *grpRecCuts = new AliDielectronCutGroup("rec","rec",AliDielectronCutGroup::kCompAND);
216   grpRecCuts->AddCut(trkRecCuts);
217   grpRecCuts->AddCut(varRecCuts);
218   die->GetTrackFilter().AddCuts(grpRecCuts);
219   grpRecCuts->Print();
220
221
222   /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv PID CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
223   AliDielectronVarCuts *pidVarCuts = new AliDielectronVarCuts("varPIDCuts","varPIDCuts");
224   pidVarCuts->AddCut(AliDielectronVarManager::kTPCnSigmaEle,  -2. ,    3.     ); //-3.0
225   pidVarCuts->AddCut(AliDielectronVarManager::kTPCnSigmaPio,   3.5, 1000.     );
226   pidVarCuts->AddCut(AliDielectronVarManager::kTPCnSigmaPro,   4. , 1000.     ); //3.0
227   //  AliDielectronPID *pidCuts        = new AliDielectronPID("PIDCuts","PIDCuts"); //not used
228   AliDielectronCutGroup *grpPIDCuts = new AliDielectronCutGroup("PID","PID",AliDielectronCutGroup::kCompAND);
229   grpPIDCuts->AddCut(pidVarCuts);
230   //grpPIDCuts->AddCut(pidCuts);
231   die->GetTrackFilter().AddCuts(grpPIDCuts);
232   grpPIDCuts->Print();
233
234   //
235
236
237   //exclude conversion electrons selected by the tender
238   //   AliDielectronTrackCuts *noconv=new AliDielectronTrackCuts("noConv","noConv");
239   //   noconv->SetV0DaughterCut(AliPID::kElectron,kTRUE);
240   //   cuts->AddCut(noconv);
241 }
242
243 //______________________________________________________________________________________
244 void SetupPairCuts(AliDielectron *die, Int_t cutDefinition)
245 {
246   //
247   // Setup the pair cuts
248   //
249
250   // add conversion rejection
251   AliDielectronVarCuts *gammaCut=new AliDielectronVarCuts("gammaCut","gammaCut");
252   gammaCut->AddCut(AliDielectronVarManager::kM,0.,.05);
253   die->GetPairPreFilter().AddCuts(gammaCut);
254
255 }
256
257 //______________________________________________________________________________________
258 void InitHistograms(AliDielectron *die, Int_t cutDefinition)
259 {
260   //
261   // Initialise the histograms
262   //
263
264   //Setup histogram Manager
265   AliDielectronHistos *histos=
266   new AliDielectronHistos(die->GetName(),
267                           die->GetTitle());
268
269   //Initialise histogram classes
270   histos->SetReservedWords("Track;Pair");
271
272   //Track classes
273   //to fill also track info from 2nd event loop until 2
274   for (Int_t i=0; i<2; ++i){
275     histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i)));
276   }
277
278   //Pair classes
279   // to fill also mixed event histograms loop until 10
280   for (Int_t i=0; i<3; ++i){
281     histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i)));
282   }
283
284   //add MC signal histograms to track and pair class
285   if(die->GetMCSignals()) {
286     for(Int_t isig=0; isig<die->GetMCSignals()->GetEntriesFast(); isig++) {
287       TString sigMCname = die->GetMCSignals()->At(isig)->GetName(); 
288
289       // mc truth
290       histos->AddClass(Form("Pair_%s_MCtruth",       sigMCname.Data()));
291       histos->AddClass(Form("Track_Legs_%s_MCtruth", sigMCname.Data())); 
292       // mc reconstructed
293       histos->AddClass(Form("Pair_%s",               sigMCname.Data()));
294       histos->AddClass(Form("Track_Legs_%s",         sigMCname.Data())); 
295     }
296   }
297
298   //add histograms to event class
299   if (cutDefinition==0) {
300     histos->AddClass("Event");
301     histos->UserHistogram("Event","","",300,-15.,15.,AliDielectronVarManager::kZvPrim);
302     histos->UserHistogram("Event","","",
303                           100,-2.,2.,100,-2.,2.,AliDielectronVarManager::kXvPrim,AliDielectronVarManager::kYvPrim);
304     histos->UserHistogram("Event","","",GetRunNumbers(),AliDielectronVarManager::kRunNumber);
305   }
306
307   //add histograms to Track classes
308   histos->UserHistogram("Track","","",200,0,20.,AliDielectronVarManager::kPt);
309   histos->UserHistogram("Track","","",
310                         400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE);
311
312   histos->UserHistogram("Track","","",
313                         100,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE);
314   histos->UserHistogram("Track","","",
315                         100,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio,kTRUE);
316   histos->UserHistogram("Track","","",
317                         100,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPro,kTRUE);
318
319   histos->UserHistogram("Track","","",
320                         100,-2,2,144,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi);
321   histos->UserHistogram("Track","","",
322                         160,-0.5,159.5,AliDielectronVarManager::kNclsTPC);
323   histos->UserHistogram("Track","","",
324                         100,0.,10.,AliDielectronVarManager::kTPCchi2Cl);
325   histos->UserHistogram("Track","","",
326                         150,-15,15,160,-0.5,159.5,AliDielectronVarManager::kYsignedIn,AliDielectronVarManager::kTPCsignalN);
327   histos->UserHistogram("Track","","",
328                         1000,0.,0.,AliDielectronVarManager::kKinkIndex0);
329
330   histos->UserHistogram("Track","","",500,-1.,1.,AliDielectronVarManager::kImpactParXY);
331   histos->UserHistogram("Track","","",600,-3.,3.,AliDielectronVarManager::kImpactParZ);
332
333   //add histograms to Pair classes
334   histos->UserHistogram("Pair","","",
335                         301,-.01,6.01,AliDielectronVarManager::kM);
336   histos->UserHistogram("Pair","","",
337                         100,-1.,1.,AliDielectronVarManager::kY);
338   histos->UserHistogram("Pair","","",
339                         100,0.,3.15,AliDielectronVarManager::kOpeningAngle);
340
341   die->SetHistogramManager(histos);
342
343 }
344
345 //______________________________________________________________________________________
346 void InitCF(AliDielectron* die, Int_t cutDefinition)
347 {
348   //
349   // Setupd the CF Manager if needed
350   //
351
352   AliDielectronCF *cf=new AliDielectronCF(die->GetName(),die->GetTitle());
353
354   //pair variables
355   cf->AddVariable(AliDielectronVarManager::kPt,"0.0, 1.0, 1.3, 2.0, 3.0, 5., 7.0, 10.0, 100.0");
356
357   cf->AddVariable(AliDielectronVarManager::kZvPrim,"-10,-5,0,5,10");
358   cf->AddVariable(AliDielectronVarManager::kY,"-1,-0.9,-0.8,-0.3,0,0.3,0.9,1.0");
359   cf->AddVariable(AliDielectronVarManager::kM,125,0.,125*.04); //40Mev Steps
360   cf->AddVariable(AliDielectronVarManager::kPairType,11,0,11);
361
362   //leg variables
363   cf->AddVariable(AliDielectronVarManager::kPt,"0.0, 0.8, 1.0, 1.1, 1.2, 1.3, 100.0",kTRUE);
364   cf->AddVariable(AliDielectronVarManager::kEta,"-1.,-0.9,-0.8,0,0.8,0.9,1.0",kTRUE);
365
366   //cf->AddVariable(AliDielectronVarManager::kITSLayerFirstCls,6,0,6,kTRUE);
367
368   if (hasMC && 0){ //ATTENTION SWITCHED OFF
369     cf->AddVariable(AliDielectronVarManager::kPdgCode,10000,-5000.5,4999.5,kTRUE);
370     cf->AddVariable(AliDielectronVarManager::kPdgCodeMother,10000,-5000.5,4999.5,kTRUE);
371     cf->AddVariable(AliDielectronVarManager::kPdgCodeGrandMother,10000,-5000.5,4999.5,kTRUE);
372   }
373
374   if(hasMC) {
375     //only steps for efficiencies
376     cf->SetStepsForMCtruthOnly();
377     //only in this case write MC truth info
378     if (cutDefinition==0){
379       cf->SetStepForMCtruth();
380     }
381   }
382   // cf->SetStepsForSignal();
383
384   die->SetCFManagerPair(cf);
385 }
386
387 //______________________________________________________________________________________
388 void SetupMCsignals(AliDielectron *die){
389   AliDielectronSignalMC* promptJpsi = new AliDielectronSignalMC("promptJpsi","Prompt J/psi");   // prompt J/psi (not from beauty decays)
390   promptJpsi->SetLegPDGs(11,-11);
391   promptJpsi->SetMotherPDGs(443,443);
392   promptJpsi->SetGrandMotherPDGs(503,503,kTRUE,kTRUE);   // not from beauty hadrons
393   promptJpsi->SetMothersRelation(AliDielectronSignalMC::kSame);
394   promptJpsi->SetFillPureMCStep(kTRUE);
395   promptJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState);
396   promptJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE);
397   promptJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE);
398   promptJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE);
399   die->AddSignalMC(promptJpsi);
400
401   AliDielectronSignalMC* promptJpsiNonRad = new AliDielectronSignalMC("promptJpsiNonRad","Prompt J/psi non-Radiative");   // prompt J/psi (not from beauty decays)
402   promptJpsiNonRad->SetLegPDGs(11,-11);
403   promptJpsiNonRad->SetMotherPDGs(443,443);
404   promptJpsiNonRad->SetGrandMotherPDGs(503,503,kTRUE,kTRUE);   // not from beauty hadrons
405   promptJpsiNonRad->SetMothersRelation(AliDielectronSignalMC::kSame);
406   promptJpsiNonRad->SetFillPureMCStep(kTRUE);
407   promptJpsiNonRad->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState);
408   promptJpsiNonRad->SetCheckBothChargesLegs(kTRUE,kTRUE);
409   promptJpsiNonRad->SetCheckBothChargesMothers(kTRUE,kTRUE);
410   promptJpsiNonRad->SetCheckBothChargesGrandMothers(kTRUE,kTRUE);
411   promptJpsiNonRad->SetJpsiRadiative(AliDielectronSignalMC::kIsNotRadiative);
412   die->AddSignalMC(promptJpsiNonRad);
413
414   AliDielectronSignalMC* promptJpsiRad = new AliDielectronSignalMC("promptJpsiRad","Prompt J/psi Radiative");   // prompt J/psi (not from beauty decays)
415   promptJpsiRad->SetLegPDGs(11,-11);
416   promptJpsiRad->SetMotherPDGs(443,443);
417   promptJpsiRad->SetGrandMotherPDGs(503,503,kTRUE,kTRUE);   // not from beauty hadrons
418   promptJpsiRad->SetMothersRelation(AliDielectronSignalMC::kSame);
419   promptJpsiRad->SetFillPureMCStep(kTRUE);
420   promptJpsiRad->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState);
421   promptJpsiRad->SetCheckBothChargesLegs(kTRUE,kTRUE);
422   promptJpsiRad->SetCheckBothChargesMothers(kTRUE,kTRUE);
423   promptJpsiRad->SetCheckBothChargesGrandMothers(kTRUE,kTRUE);
424   promptJpsiRad->SetJpsiRadiative(AliDielectronSignalMC::kIsRadiative);
425   die->AddSignalMC(promptJpsiRad);
426
427 }
428
429 TVectorD *GetRunNumbers() {
430   // returns a vector with the runnumber used in the period
431
432   Double_t first=0;
433   Double_t last =1;
434
435   switch(iPeriod) {
436   case k10b: first=114737; last=117223; break;
437   case k10c: first=117777; last=121417; break;
438   case k10d: first=121692; last=126437; break;
439   case k10e: first=127102; last=130850; break;
440   case k10f: first=130931; last=135031; break;
441   case k10h: first=136831; last=139517; break;
442   case k11a: first=141052; last=146974; break;
443   case k11d: first=155838; last=159649; break;
444   case k11h: first=165772; last=170718; break;
445   case k12h: first=188720; last=192738; break;
446   }
447   //  printf("iPeriod: %d \t %.0f-%.0f \n",iPeriod,first,last);
448   return (AliDielectronHelper::MakeLinBinning(last-first, first, last));
449 }