Analysis task to study the extra SPD clusters.
[u/mrichter/AliRoot.git] / PWGPP / forward / SPDClustTask / AnalysisSPDClustTask.C
1 Bool_t needRecPoints = kTRUE;
2
3 void AnalysisSPDClustTask
4 (
5  TString dataset="/alice/sim/LHC11d2_000119161",
6  TString outFName = "trbg.root",
7  Int_t   nEvents   = -1,
8  Float_t etaMin     =-1.5,        // min eta range to fill in histos
9  Float_t etaMax     = 1.5,        // max eta range to fill in histos
10  Float_t zMin       = -5,         // process events with Z vertex min
11  Float_t zMax       =  5,         //                     max positions
12  Float_t cutSigNStd  = 1.,        // cut on weighed distance used to extract signal
13  Float_t cutSigDPhiS = -1,        // cut on dPhi-phiBent used to extract signal (if negative -> dphi*sqrt(cutSigNStd)
14  Bool_t  useMC  = kTRUE,          // fill MC info (doRec=kTRUE)
15  // specific parameters for reconstruction
16  Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
17  float  nStdDev     = 1.,         // number of st.dev. for tracklet cut to keep
18  float  dphi        = 0.06,        // dphi window (sigma of tracklet cut)
19  float  dtht        = 0.025,       // dtheta .... (if negative, abs will be used with additional cut on |dthetaX|, apart from w.distance
20  float  phishift    = 0.0045,      // bending shift
21  Bool_t remOvl      = kTRUE,       
22  float  ovlPhiCut   = 0.005, 
23  float  ovlZetaCut  = 0.05,
24  Int_t  nEventsSkip = 0
25  )
26 {
27   //  
28   if (cutSigDPhiS<0) cutSigDPhiS = TMath::Sqrt(cutSigNStd)*dphi;
29   //
30   printf("Start Analysis for %s, max %d Events skipping %d, Event Cuts: %.1f<eta<%.1f, %.2f<Zv<%.2f\n",
31          dataset.Data(),nEvents,nEventsSkip,etaMin,etaMax,zMin,zMax);
32   printf("Tracklet cuts: dPhi:%.3f dTheta:%.3f phiShift:%.4f | Keep %.1f NstDev\n"
33          "Scale dTheta: %s | Signal Selection: NstDev:%.1f, dPhiS: %.3f\n", 
34          dphi,dtht,phishift,nStdDev,scaleDTheta ? "ON":"OFF",
35          cutSigNStd,cutSigDPhiS);
36   //
37   if (nEvents<0) nEvents = int(1e9);
38   TString format = GetFormatFromDataSet(dataset);
39   //
40   // ALICE stuff
41   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
42   if (!mgr) mgr = new AliAnalysisManager("Test train");
43   //
44   InputHandlerSetup(format,useMC);
45   // compile our task
46   gProof->Load("AliSPDClustTask.cxx++");
47   //
48   // load and run AddTask macro
49   //  gROOT->LoadMacro("AddMultTaskTrackletMulti.C");
50   //
51   // create our task
52   
53   AliSPDClustTask *task = new AliSPDClustTask("AliSPDClustTask");
54   // create output container
55   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist", TList::Class(),AliAnalysisManager::kOutputContainer,outFName);
56   //  coutput1->SetSpecialOutput();
57   // add our task to the manager
58   mgr->AddTask(task);
59
60   // finaly connect input and output
61   mgr->ConnectInput(task, 0,  mgr->GetCommonInputContainer());
62   mgr->ConnectOutput(task,1,coutput1);
63   //  mgr->SetSpecialOutputLocation("root://alicers01.cern.ch//tmp/output/");
64   //
65   //
66   task->SetUseMC(useMC);
67   //
68   task->SetEtaMin(etaMin);
69   task->SetEtaMax(etaMax);
70   task->SetZVertexMin(zMin);
71   task->SetZVertexMax(zMax);
72   //
73   task->SetDPhiSCut(cutSigDPhiS);
74   task->SetNStdCut(cutSigNStd);
75   //
76   task->SetScaleDThetaBySin2T(scaleDTheta);
77   task->SetNStdDev(nStdDev);
78   task->SetPhiWindow(dphi);
79   task->SetThetaWindow(dtht);
80   task->SetPhiShift(phishift);
81   task->SetPhiOverlapCut(ovlPhiCut);
82   task->SetZetaOverlapCut(ovlZetaCut);
83   task->SetRemoveOverlaps(remOvl);
84   //
85   task->SetInput("spectraCombine.root");
86   //
87   printf("Requesting physics selection in %s mode\n",useMC ? "MC":"Data");
88   gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
89   //  /*
90   //gROOT->ProcessLine(".L AddTaskPhysicsSelection.C");
91   AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(useMC);
92   task->SelectCollisionCandidates();//AliVEvent::kMB);
93   //
94   //  */
95   // Run analysis
96   mgr->InitAnalysis();
97   // process dataset  
98   mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip); 
99   //
100   TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_"; 
101   evstCmd += outFName;  evstCmd += " \nfi";
102   gSystem->Exec( evstCmd.Data() );
103   
104 }
105
106
107 TString GetFormatFromDataSet(TString dataset) {
108   
109 //   Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
110   TString dsTreeName;
111   if (dataset.Contains("#")) {
112     Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
113     dsTreeName=dataset(dataset.Last('#'),dataset.Length());
114   } else {
115     Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
116     TFileCollection *ds = gProof->GetDataSet(dataset.Data());
117     if (!ds) {
118       Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
119       return "";
120     }
121     dsTreeName = ds->GetDefaultTreeName();
122   }
123
124   if (dsTreeName.Contains("esdTree")) {
125     Info("runAAF.C","ESD input format detected ...");
126     return "ESD";
127   } else if (dsTreeName.Contains("aodTree"))  {
128     Info("runAAF.C","AOD input format detected ...");
129     return "AOD";
130   } else {
131     Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
132     Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
133   }
134   
135   return "";
136 }
137
138 Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
139 {
140   format.ToLower();
141
142   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
143
144   AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
145
146   if (cin) return;
147
148   if (!format.CompareTo("esd"))
149   {
150     AliESDInputHandler *esdInputHandler = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
151
152     if (!esdInputHandler)
153     {
154       Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
155       if (needRecPoints)
156         esdInputHandler = new AliESDInputHandlerRP();
157       else 
158         esdInputHandler = new AliESDInputHandler();
159       //
160       mgr->SetInputEventHandler(esdInputHandler);
161     }
162     //
163     if (useKine)
164     {
165       AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
166
167       if (!mcInputHandler)
168       {
169         Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
170         AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
171         mgr->SetMCtruthEventHandler(mcInputHandler);
172       }
173     }
174
175   }
176   else if (!format.CompareTo("aod"))
177   {
178     AliAODInputHandler *aodInputHandler = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
179
180     if (!aodInputHandler)
181     {
182       Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
183       aodInputHandler = new AliAODInputHandler();
184       mgr->SetInputEventHandler(aodInputHandler);
185     }
186   }
187   else
188   {
189     Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
190     return kFALSE;
191   }
192
193   return kTRUE;
194 }
195
196 void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
197                      float zMin, float zMax, float zMixBinSz)
198 {
199   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
200   if (!mgr) return;
201   int bufferSize = 1;
202   AliESDInputHandlerRP *esdH = dynamic_cast<AliESDInputHandlerRP*>(mgr->GetInputEventHandler());
203   if (!esdH) return;
204   //
205   AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
206   esdMixH->SetInputHandlerForMixing(esdH);
207   AliMixEventPool *evPool = new AliMixEventPool("Pool");
208   AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
209   AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
210   //  evPool->AddCut(tracklets);
211   evPool->AddCut(zvertex);
212   //evPool->Init();
213   evPool->Print();
214   esdMixH->SetEventPool(evPool);
215   esdH->SetMixingHandler(esdMixH);
216 }
217
218 void AddPhysicsSelection(Bool_t isMC)
219 {
220   // physics selection a la Michele
221   if(!isMC) {
222     //AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
223     //    physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
224     /*
225     physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
226     physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
227     */
228     //
229     //    physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
230     //    physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
231     //
232     // This are needed only to fill the statistics tables
233     //    physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
234     //    physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
235     //    physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
236     //
237     /*
238     physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
239     physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
240     physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
241     //
242     physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
243     physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
244     physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
245     //
246     */
247     /*
248     //
249     physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
250     physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
251     physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
252     //
253     physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
254     physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
255     physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
256     */
257   } 
258   // if you use the following line, your task only gets the selected events
259   //  task->SelectCollisionCandidates(AliVEvent::kUserDefined);
260   //  task->SelectCollisionCandidates();
261   //
262   //Alternatively, in the UserExec of your task:
263   //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
264   //
265 }