/[PAMELA software]/PamelaLevel2/src/PamLevel2.cpp
ViewVC logotype

Contents of /PamelaLevel2/src/PamLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.34 - (show annotations) (download)
Wed Mar 21 10:49:43 2007 UTC (17 years, 8 months ago) by pam-fi
Branch: MAIN
Changes since 1.33: +16 -9 lines
...

1 #include <PamLevel2.h>
2 //--------------------------------------
3 //
4 //
5 //--------------------------------------
6 /**
7 * Default constructor
8 */
9 PamTrack::PamTrack(){
10 trk_track = 0;
11 calo_track = 0;
12 tof_track = 0;
13 candeleteobj = 0;
14 };
15 //--------------------------------------
16 //
17 //
18 //--------------------------------------
19 /**
20 * Constructor
21 */
22 PamTrack::PamTrack(TrkTrack* t, CaloTrkVar* c, ToFTrkVar* o){
23
24 trk_track = 0;
25 calo_track = 0;
26 tof_track = 0;
27 // if(t)trk_track = new TrkTrack(*t);
28 // if(c)calo_track = new CaloTrkVar(*c);
29 // if(o)tof_track = new ToFTrkVar(*o);
30 if(t)trk_track = t;
31 if(c)calo_track = c;
32 if(o)tof_track = o;
33 candeleteobj = 0;
34
35 };
36 PamTrack::PamTrack(const PamTrack& track){
37
38 TrkTrack *t = track.trk_track;
39 CaloTrkVar *c = track.calo_track;
40 ToFTrkVar *o = track.tof_track;
41
42 trk_track = 0;
43 calo_track = 0;
44 tof_track = 0;
45 if(t)trk_track = new TrkTrack(*t);
46 if(c)calo_track = new CaloTrkVar(*c);
47 if(o)tof_track = new ToFTrkVar(*o);
48 candeleteobj = 1;
49
50 }
51 void PamTrack::Clear(){
52
53 // cout << "PamTrack::Clear() "<<candeleteobj<<endl;
54 if(candeleteobj){
55 if(trk_track) trk_track->TrkTrack::Clear();
56 if(calo_track) calo_track->CaloTrkVar::Clear();//???
57 if(tof_track) tof_track->ToFTrkVar::Clear();//???
58 }else{
59 trk_track = 0;
60 calo_track = 0;
61 tof_track = 0;
62 }
63 }
64 void PamTrack::Delete(){
65 // cout << "PamTrack::Delete() "<<candeleteobj<<endl;
66 if(candeleteobj){
67 if(trk_track) {
68 trk_track->TrkTrack::Clear();
69 delete trk_track;
70 }
71 if(calo_track){
72 calo_track->CaloTrkVar::Clear();//???
73 delete calo_track;
74 }
75 if(tof_track){
76 tof_track->ToFTrkVar::Clear();//???
77 delete tof_track;
78 }
79 }else{
80 Clear();
81 }
82 }
83 //--------------------------------------
84 //
85 //
86 //--------------------------------------
87 /**
88 * Constructor
89 */
90 PamLevel2::PamLevel2(){
91
92 // trk2_obj = TrkLevel2::GetTrkLevel2();
93 // trk1_obj = TrkLevel1::GetTrkLevel1();
94 // trkh_obj = TrkHough::GetTrkHough();
95 // calo1_obj = CaloLevel1::GetCaloLevel1();
96 // calo2_obj = CaloLevel2::GetCaloLevel2();
97 // tof_obj = ToFLevel2::GetToFLevel2();
98 // trig_obj = TrigLevel2::GetTrigLevel2();
99 // s4_obj = S4Level2::GetS4Level2();
100 // nd_obj = NDLevel2::GetNDLevel2();
101 // ac_obj = AcLevel2::GetAcLevel2();
102 // orb_obj = OrbitalInfo::GetOrbitalInfo();
103
104 h0_obj = 0;
105 trk0_obj = 0;
106 calo0_obj = 0;
107
108 trk2_obj = 0;
109 trk1_obj = 0;
110 trkh_obj = 0;
111 calo1_obj = 0;
112 calo2_obj = 0;
113 tof_obj = 0;
114 trig_obj = 0;
115 s4_obj = 0;
116 nd_obj = 0;
117 ac_obj = 0;
118 orb_obj = 0;
119
120 run_obj = 0;//new GL_RUN();
121 soft_obj = 0;// Emiliano
122 irun = -1LL;
123 runfirstentry = 0LL;
124 runlastentry = 0LL;
125
126 l0_file = NULL;
127 l0_tree = NULL;
128 iroot = -1;
129 dbc = 0;
130
131 prevshift = 0;
132
133 run_tree = NULL;
134 run_tree_clone = NULL;
135 sel_tree = NULL;
136 sel_tree_clone = NULL;
137
138 irunentry = -1LL;
139 pam_tree = NULL;
140 for(Int_t i=0; i<NCLONES; i++ )pam_tree_clone[i]=NULL;
141
142 host = "mysql://localhost/pamelaprod";
143 user = "anonymous";
144 psw = "";
145 const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
146 const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
147 const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
148 if ( !pamdbhost ) pamdbhost = "";
149 if ( !pamdbuser ) pamdbuser = "";
150 if ( !pamdbpsw ) pamdbpsw = "";
151 if ( strcmp(pamdbhost,"") ) host = pamdbhost;
152 if ( strcmp(pamdbuser,"") ) user = pamdbuser;
153 if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
154
155
156
157 // sorted_tracks = 0;//new TRefArray();
158
159 CAL0 = false;
160 CAL1 = true;
161 CAL2 = true;
162 TRK2 = true;
163 TRK1 = false;
164 TRKh = false;
165 TRKh = false;
166 TRG = true;
167 TOF = true;
168 TOF0 = false;
169 S4 = true;
170 ND = true;
171 AC = true;
172 ORB = true;
173
174 RUN = true;
175
176 SELLI = -1;
177
178 tsorted=0;
179 timage=0;
180
181 };
182 /**
183 * Delete the event (act as Dtor)
184 */
185 void PamLevel2::Delete(){
186
187 if(run_obj)delete run_obj;
188 if(soft_obj)delete soft_obj; //Emiliano
189
190 // cout << "void PamLevel2::Clear()"<<endl;
191 if(h0_obj) delete h0_obj;
192 if(trk0_obj) delete trk0_obj;
193 if(calo0_obj) delete calo0_obj;
194 if(trk1_obj) delete trk1_obj;
195 if(trk2_obj) delete trk2_obj;
196 if(trkh_obj) delete trkh_obj;
197 if(calo1_obj)delete calo1_obj;
198 if(calo2_obj)delete calo2_obj;
199 if(tof_obj) delete tof_obj;
200 if(trig_obj) delete trig_obj;
201 if(s4_obj) delete s4_obj;
202 if(nd_obj) delete nd_obj;
203 if(ac_obj) delete ac_obj;
204 if(orb_obj) delete orb_obj;
205
206 if(tsorted){
207 tsorted->Delete();
208 delete tsorted;
209 }
210 if(timage){
211 timage->Delete();
212 delete timage;
213 }
214
215 if(dbc){
216 dbc->Close();
217 delete dbc;
218 }
219
220 if(l0_file)l0_file->Close();
221 if(pam_tree)pam_tree->Delete();;
222 if(run_tree)run_tree->Delete();;
223 if(sel_tree)sel_tree->Delete();;
224 for(Int_t i=0; i<NCLONES; i++ )if(pam_tree_clone[i])pam_tree_clone[i]->Delete();;
225 if(run_tree_clone)run_tree_clone->Delete();;
226 if(sel_tree_clone)sel_tree_clone->Delete();;
227
228
229 };
230
231 /**
232 * Clear the event (NB! does not deallocate objects)
233 */
234 void PamLevel2::Clear(){
235
236 // cout << "void PamLevel2::Clear()"<<endl;
237
238 //
239 // This method is called once for every entry but RunInfo and SoftInfo do not change until the next run so we cannot clear them here unless we don't
240 // want to load them for each event even if they are the same...
241 //
242 // if(run_obj)delete run_obj;
243 // if(run_obj) run_obj->Clear(); // Emiliano: Do not deallocate run_obj here, it will give segmentation fault! call clear instead
244 // if(soft_obj) soft_obj->Clear();
245
246 if(h0_obj) h0_obj->Clear();
247 // if(trk0_obj) trk0_obj->Clear();
248 if(calo0_obj) calo0_obj->Clear();
249 if(trk1_obj) trk1_obj->Clear();
250 if(trk2_obj) trk2_obj->Clear();
251 if(trkh_obj) trkh_obj->Clear();
252 if(calo1_obj)calo1_obj->Clear();
253 if(calo2_obj)calo2_obj->Clear();
254 if(tof_obj) tof_obj->Clear();
255 if(trig_obj) trig_obj->Clear();
256 if(s4_obj) s4_obj->Clear();
257 if(nd_obj) nd_obj->Clear();
258 if(ac_obj) ac_obj->Clear();
259 if(orb_obj) orb_obj->Clear();
260
261 // if(sorted_tracks)sorted_tracks->Clear();
262 // sorted_tracks.Clear();
263
264 if(tsorted){
265 tsorted->Delete();
266 }
267 if(timage){
268 timage->Delete();
269 }
270
271 };
272
273
274 //--------------------------------------
275 //
276 //
277 //--------------------------------------
278 void *PamLevel2::GetPointerTo(const char* c ){
279
280 TString objname = c;
281
282 if(!objname.CompareTo("TrkLevel1")) {
283 if(!trk1_obj){
284 trk1_obj = new TrkLevel1();
285 trk1_obj->Set();
286 }
287 return &trk1_obj;
288 };
289 if(!objname.CompareTo("TrkLevel2")) {
290 if(!trk2_obj){
291 trk2_obj = new TrkLevel2();
292 trk2_obj->Set();
293 }
294 return &trk2_obj;
295 };
296 if(!objname.CompareTo("TrkHough")) {
297 if(!trkh_obj) trkh_obj = new TrkHough();
298 return &trkh_obj;
299 };
300 if(!objname.CompareTo("CaloLevel1")) {
301 if(!calo1_obj) calo1_obj = new CaloLevel1();
302 return &calo1_obj;
303 };
304 if(!objname.CompareTo("CaloLevel2")) {
305 if(!calo2_obj){
306 calo2_obj = new CaloLevel2();
307 calo2_obj->Set();
308 };
309 return &calo2_obj;
310 };
311 if(!objname.CompareTo("ToFLevel2")) {
312 if(!tof_obj){
313 tof_obj = new ToFLevel2();
314 tof_obj->Set();
315 }
316 return &tof_obj;
317 };
318 if(!objname.CompareTo("TrigLevel2")) {
319 if(!trig_obj) trig_obj = new TrigLevel2();
320 return &trig_obj;
321 };
322 if(!objname.CompareTo("S4Level2")) {
323 if(!s4_obj) s4_obj = new S4Level2();
324 return &s4_obj;
325 };
326 if(!objname.CompareTo("NDLevel2")) {
327 if(!nd_obj) nd_obj = new NDLevel2();
328 return &nd_obj;
329 };
330 if(!objname.CompareTo("AcLevel2")) {
331 if(!ac_obj) ac_obj = new AcLevel2();
332 return &ac_obj;
333 };
334 if(!objname.CompareTo("OrbitalInfo")){
335 if(!orb_obj) orb_obj = new OrbitalInfo();
336 return &orb_obj;
337 };
338
339 if(!objname.CompareTo("RunInfo"))return &run_obj;
340
341 if(!objname.CompareTo("SoftInfo"))return &soft_obj; // Emiliano
342
343 return NULL;
344 };
345 //--------------------------------------
346 //
347 //
348 //--------------------------------------
349 /**
350 * Retrieves the calorimeter track matching the seqno-th tracker stored track.
351 * (If seqno = -1 retrieves the self-trigger calorimeter track)
352 */
353 CaloTrkVar *PamLevel2::GetCaloStoredTrack(int seqno){
354
355 if( !calo2_obj )return 0;
356
357 if( calo2_obj->CaloLevel2::ntrk()==0 ){
358 cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no Calorimeter tracks are stored"<<endl;
359 return NULL;
360 };
361
362 CaloTrkVar *c = 0;
363 Int_t it_calo=0;
364
365 do{
366 c = calo2_obj->CaloLevel2::GetCaloTrkVar(it_calo);
367 it_calo++;
368 } while( c && seqno != c->trkseqno && it_calo < calo2_obj->CaloLevel2::ntrk());
369
370 if(!c || seqno != c->trkseqno){
371 c = 0;
372 if(seqno!=-1)cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match Calorimeter stored tracks"<<endl;
373 };
374 return c;
375
376 };
377 //--------------------------------------
378 //
379 //
380 //--------------------------------------
381 /**
382 * Retrieves the ToF track matching the seqno-th tracker stored track.
383 * (If seqno = -1 retrieves the tracker-independent tof track)
384 */
385 ToFTrkVar *PamLevel2::GetToFStoredTrack(int seqno){
386
387 if( !tof_obj )return 0;
388
389 if( tof_obj->ToFLevel2::ntrk()==0 ){
390 cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no ToF tracks are stored"<<endl;
391 return NULL;
392 };
393
394 ToFTrkVar *c = 0;
395 Int_t it_tof=0;
396
397 do{
398 c = tof_obj->ToFLevel2::GetToFTrkVar(it_tof);
399 it_tof++;
400 } while( c && seqno != c->trkseqno && it_tof < tof_obj->ToFLevel2::ntrk());
401
402 if(!c || seqno != c->trkseqno){
403 c = 0;
404 if(seqno!=-1)cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match ToF stored tracks"<<endl;
405 };
406 return c;
407
408 };
409
410 //--------------------------------------
411 //
412 //
413 //--------------------------------------
414 /**
415 * Give the pamela track associated to a tracker track, retrieving related calorimeter and tof track information.
416 */
417 PamTrack* PamLevel2::GetPamTrackAlong(TrkTrack* t){
418
419 cout <<"PamLevel2::GetPamTrackAlong(TrkTrack* t) **obsolete** "<<endl;
420 cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
421
422 CaloTrkVar *c = 0;
423 ToFTrkVar *o = 0;
424
425 if(CAL2) c = GetCaloStoredTrack(t->GetSeqNo());
426 if(TOF) o = GetToFStoredTrack(t->GetSeqNo());
427
428 // if(t && c && o)track = new PamTrack(t,c,o);
429 PamTrack *track = new PamTrack(t,c,o);
430
431 return track;
432
433 };
434 //--------------------------------------
435 //
436 //
437 //--------------------------------------
438 /**
439 * Retrieves the it-th stored track.
440 * It override TrkLevel2::GetTrack(int it).
441 * @param itrk Track number, ranging from 0 to GetNTracks().
442 */
443
444 PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){
445
446 cout <<"PamLevel2::GetStoredTrack(Int_t itrk) **to-be-updated** "<<endl;
447 cout <<"for the moment, better use separately the methods: TrkLevel2::GetStoredTrack(seqno) CaloLevel2::GetCaloTrkVar(Int_t notrack) ToFLevel2::GetToFTrkVar(Int_t notrack)"<<endl;
448 cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
449 PamTrack *track = 0;
450
451 if( itrk >=0 && itrk < trk2_obj->TrkLevel2::ntrk() ){
452
453 TrkTrack *t = trk2_obj->TrkLevel2::GetStoredTrack(itrk);
454 track = GetPamTrackAlong(t);
455
456
457
458 }else{
459 cout << "PamLevel2::GetStoredTrack(int) : tracker track SeqNo "<< itrk <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
460 };
461
462 return track;
463
464 }
465 //--------------------------------------
466 //
467 //
468 //--------------------------------------
469 /**
470 * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is PamLevel2::sorted_tracks.
471 * @param how String to set the sorting cryterium (es: "CAL" or "TRK+CAL+TOF" ecc...).
472 * Sorting cryteria:
473 * TRK: lower chi**2
474 * CAL: lower Y spatial residual on the first calorimeter plane
475 * TOF: bigger numebr of hit PMTs along the track, on S12 S21 S32 (where paddles are along the Y axis).
476 * The default sorting cryterium is "TOF+CAL".
477 *
478 * The total number of physical tracks is always given by GetNTracks() and the it-th physical track can be retrieved by means of the methods GetTrack(int it) and GetTrack(int it, TString how).
479 */
480 void PamLevel2::SortTracks(TString how){
481
482 // cout <<" PamLevel2::SortTracks(TString how) "<<endl;
483 if( !trk2_obj ){
484 cout << "void PamLevel2::SortTracks(TString how): TrkLevel2 not loaded !!!";
485 return;
486 };
487 // cout << "call SortTracs() "<<endl;
488 //Save current Object count
489 Int_t ObjectNumber = TProcessID::GetObjectCount();
490
491 // cout << "ObjectNumber "<<ObjectNumber <<endl;
492
493 // if(!sorted_tracks)sorted_tracks = new TRefArray();
494 // sorted_tracks->Clear();
495 // sorted_tracks.Clear();
496
497 if(!tsorted)tsorted = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
498 tsorted->Delete();
499 TClonesArray &ttsorted = *tsorted;
500 if(!timage)timage = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
501 timage->Delete();
502 TClonesArray &ttimage = *timage;
503
504
505 // loop over the tracks sorted by the tracker
506 Bool_t use_TRK = how.Contains("TRK", TString::kIgnoreCase);
507 Bool_t use_CAL = how.Contains("CAL", TString::kIgnoreCase);
508 Bool_t use_TOF = how.Contains("TOF", TString::kIgnoreCase);
509
510 if( !CAL2 && use_CAL) use_CAL = false;
511 if( !TOF && use_TOF) use_TOF = false;
512
513 if( !TRK2 ){
514 // cout << "SortTracks() : without tracker does not work!!! (not yet)" << endl;
515 return;
516 };
517
518 // cout << "use_CAL "<<use_CAL<<" use_TOF "<<use_TOF<<" use_TRK "<<use_TRK <<endl;
519
520 for(Int_t i=0; i < trk2_obj->TrkLevel2::GetNTracks(); i++){
521
522 TrkTrack *ts = 0;
523 CaloTrkVar *cs = 0;
524 ToFTrkVar *os = 0;
525
526 // get tracker tracks
527 TrkTrack *tp = trk2_obj->TrkLevel2::GetTrack(i); //tracker
528 CaloTrkVar *cp = GetCaloStoredTrack(tp->GetSeqNo());
529 ToFTrkVar *op = GetToFStoredTrack(tp->GetSeqNo());
530
531 TrkTrack *ti = 0; //tracker (image)
532 CaloTrkVar *ci = 0;
533 ToFTrkVar *oi = 0;
534 // cout << "trk track n. "<<i << " "<<hex<< tp <<dec<< endl;
535 // if track has an image, check image selection
536 if(tp->HasImage()){
537
538 ti = trk2_obj->TrkLevel2::GetTrackImage(i); //tracker (image)
539 ci = GetCaloStoredTrack(ti->GetSeqNo());
540 oi = GetToFStoredTrack(ti->GetSeqNo());
541
542 // cout << "its image "<<i << " "<<hex<< ti <<dec<< endl;
543
544 //assign starting scores
545 Int_t tp_score = 0; //main track sorted by the tracker
546 Int_t ti_score = 0; //image track
547
548 // ------------------------
549 // calorimeter check
550 // ------------------------
551 // check the Y spatial residual on the first calorimeter plane
552 // (cut on calorimeter variables from Emiliano)
553 if( use_CAL && !calo2_obj ){
554 cout << "void PamLevel2::SortTracks(TString how): how= "<<how<<" but CaloLevel2 not loaded !!!";
555 return;
556 };
557 if(
558 use_CAL &&
559 calo2_obj->npcfit[1] > 15 && //no. of fit planes on Y view
560 calo2_obj->varcfit[1] < 1000. && //fit variance on Y view
561 cp && ci &&
562 true){
563
564
565 Float_t resy_p = cp->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_p < 0)resy_p= - resy_p;
566 Float_t resy_i = ci->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_i < 0)resy_i= - resy_i;
567
568 if(resy_p <= resy_i) tp_score++;
569 else ti_score++;
570
571 // cout << "CALO "<<tp_score<<ti_score<<endl;
572
573 };
574 // ------------------------
575 // TOF check
576 // ------------------------
577 // check the number of hit pmts along the track
578 // on S12 S21 and S32, where paddles are parallel to Y axis
579 if( use_TOF && !tof_obj ){
580 cout << "void PamLevel2::SortTracks(TString how): how= "<<how<<" but ToFLevel2 not loaded !!!";
581 return;
582 };
583 if( use_TOF && op && oi ){
584
585 Int_t nphit_p =0;
586 Int_t nphit_i =0;
587
588
589 /* cout << "track: npmtadc "<< op->npmtadc << endl;
590 cout << "track: npmttdc "<< op->npmttdc << endl;
591 cout << "image: npmtadc "<< oi->npmtadc << endl;
592 cout << "image: npmttdc "<< oi->npmttdc << endl;*/
593
594 for (Int_t ih=0; ih < op->npmtadc; ih++){
595 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
596 if(pl == 1 || pl == 2 || pl == 5)nphit_p++;
597 };
598
599 for (Int_t ih=0; ih < oi->npmtadc; ih++){
600 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
601 if(pl == 1 || pl == 2 || pl == 5)nphit_i++;
602 };
603
604 if(
605 use_TOF &&
606 (nphit_p+nphit_i) !=0 &&
607 true){
608
609 if( nphit_p >= nphit_i) tp_score++;
610 else ti_score++;
611 };
612 // cout << "TOF "<<tp_score<<ti_score<<endl;
613 };
614 if(tp_score == ti_score) use_TRK = true;
615 // ------------------------
616 // tracker check
617 // ------------------------
618 // chi**2 difference is not always large enough to distinguish among
619 // the real track and its image.
620 // Tracker check will be applied always when calorimeter and tof information is ambiguous.
621 if(use_TRK){
622 if( tp->chi2 > 0 && tp->chi2 < ti->chi2 ) tp_score++ ;
623 else if( ti->chi2 > 0 && ti->chi2 < tp->chi2 ) ti_score++ ;
624 // cout << "TRK "<<tp_score<<ti_score<<endl;
625 };
626
627 // ------------------------
628 // the winner is....
629 // ------------------------
630 if (tp_score > ti_score) {
631 // ts = tp;//the track sorted by the tracker!!
632 // cs = cp;
633 // os = op;
634 }else if (tp_score < ti_score) {
635 ts = ti;//its image!!
636 cs = ci;
637 os = oi;
638
639 ti = tp;//its image!!
640 ci = cp;
641 oi = op;
642
643 tp = ts;//its image!!
644 cp = cs;
645 op = os;
646
647
648 }else {
649 // ts = tp;
650 // cs = cp;
651 // os = op;
652 // cout << "Warning - track image ambiguity not solved" << endl;
653 // cout << ts->GetNtot() << " "<< ts->chi2 << " " << npcfit[1] << " "<< nphit_p << endl;
654 };
655
656 }else{
657 // ts = tp;
658 // cs = cp;
659 // os = op;
660 };
661
662 // cout <<" SortTracks() "<<i<<" -- "<<ts<<endl;
663 // sorted_tracks->Add(ts);//save the track in the sorted array
664 // sorted_tracks.Add(ts);//save the track in the sorted array
665 // sorted_tracks.Add(tp);//save the track in the sorted array
666 // cout << "SortTracks:: sorted_tracks->Add(it) "<<i<<" "<<ts<<endl;
667 // cout<<"o "<<tp<<endl;
668 // cout<<"o "<<cp<<endl;
669 // cout<<"o "<<op<<endl;
670 new(ttsorted[i]) PamTrack(tp,cp,op);
671 new(ttimage[i]) PamTrack(ti,ci,oi);
672 };
673
674 if( tsorted->GetEntries() != trk2_obj->GetNTracks() ){
675 cout << "void PamLevel2::SortTracks(TString how): tsorted->GetEntries() "<<tsorted->GetEntries()<<" != trk2_obj->GetNTracks() = "<<trk2_obj->GetNTracks() <<endl;
676 tsorted->Delete(); tsorted=0;
677 timage->Delete(); timage=0;
678 }
679
680 //Restore Object count
681 //To save space in the table keeping track of all referenced objects
682 //We reset the object count to what it was at the beginning of the event.
683 TProcessID::SetObjectCount(ObjectNumber);
684
685 };
686 //--------------------------------------
687 //
688 //
689 //--------------------------------------
690 /**
691 * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.
692 * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.
693 */
694 // TRefArray *PamLevel2::GetTracks(){
695
696 // // *-*-*-*-*-*-*-*-*-*-*-*-*
697 // SortTracks("+CAL+TOF");
698 // // *-*-*-*-*-*-*-*-*-*-*-*-*
699
700 // // return sorted_tracks;
701 // return &sorted_tracks;
702
703 // };
704 TClonesArray *PamLevel2::GetTracks(){
705
706 // *-*-*-*-*-*-*-*-*-*-*-*-*
707 SortTracks("+CAL+TOF");
708 // *-*-*-*-*-*-*-*-*-*-*-*-*
709
710 return tsorted;
711
712 };
713 //--------------------------------------
714 //
715 //
716 //--------------------------------------
717 /**
718 * Retrieves the it-th Pamela "physical" track.
719 * It override TrkLevel2::GetTrack(int it).
720 * @param it Track number, ranging from 0 to GetNTracks().
721 */
722 PamTrack *PamLevel2::GetTrack(int it){
723
724 // if(!trk2_obj) return 0;
725
726 // // *-*-*-*-*-*-*-*-*-*-*-*-*
727 // SortTracks("+CAL+TOF");
728 // // *-*-*-*-*-*-*-*-*-*-*-*-*
729 // // if(!sorted_tracks)return 0;
730 // if(sorted_tracks.IsEmpty())return 0;
731
732 // PamTrack *track = 0;
733
734 // if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() && it<sorted_tracks.GetEntriesFast() ){
735 // TrkTrack *t = (TrkTrack*)sorted_tracks.At(it);
736 // track = GetPamTrackAlong(t);
737 // }else{
738 // cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
739 // };
740
741 // return track;
742
743 // cout << "PamLevel2::GetTrack(int it) "<<endl;
744 // *-*-*-*-*-*-*-*-*-*-*-*-*
745 SortTracks("+CAL+TOF");
746 // *-*-*-*-*-*-*-*-*-*-*-*-*
747 if(!tsorted)return 0;
748 if(!tsorted->GetEntries())return 0;
749
750 PamTrack *track = 0;
751
752 if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks()){
753 // TrkTrack *t = (TrkTrack*)sorted_tracks.At(it);
754 // track = GetPamTrackAlong(t);
755 TClonesArray &t = *(tsorted);
756 track = (PamTrack*)t[it];
757
758 }else{
759 cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
760 };
761
762 return track;
763
764 };
765
766 //--------------------------------------
767 //
768 //
769 //--------------------------------------
770 /**
771 * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::SortTracks().
772 * @param it Track number, ranging from 0 to GetNTracks().
773 */
774 PamTrack *PamLevel2::GetTrackImage(int it){
775
776 // if(!trk2_obj) return 0;
777
778 // // *-*-*-*-*-*-*-*-*-*-*-*-*
779 // SortTracks("+CAL+TOF");
780 // // *-*-*-*-*-*-*-*-*-*-*-*-*
781 // // if(!sorted_tracks)return 0;
782 // if(sorted_tracks.IsEmpty())return 0;
783
784 // PamTrack *image = 0;
785
786 // if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
787 // TrkTrack *temp = (TrkTrack*)sorted_tracks.At(it);
788 // if( temp->HasImage() ){
789 // TrkTrack *t = trk2_obj->TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());
790 // image = GetPamTrackAlong(t);
791 // }else{
792 // cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
793 // };
794 // }else{
795 // cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
796 // };
797
798 // return image;
799
800
801 // *-*-*-*-*-*-*-*-*-*-*-*-*
802 SortTracks("+CAL+TOF");
803 // *-*-*-*-*-*-*-*-*-*-*-*-*
804 if(!timage)return 0;
805 if(!timage->GetEntries())return 0;
806
807 PamTrack *image = 0;
808
809 if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
810 TClonesArray &t = *(tsorted);
811 PamTrack *temp = (PamTrack*)t[it];
812 if( temp->GetTrkTrack()->HasImage() ){
813 TClonesArray &t = *(timage);
814 image = (PamTrack*)t[it];
815 }else{
816 // cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
817 };
818 }else{
819 cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
820 };
821
822 return image;
823 }
824
825 //--------------------------------------
826 //
827 //
828 //--------------------------------------
829 /**
830 * Get the Pamela detector trees in a single file and make them friends.
831 * @param f TFile pointer
832 * @param detlist String to select trees to be included
833 * @return Pointer to a TTree
834 */
835 TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ){
836
837 // if( !detlist.IsNull() ) SetWhichTrees(detlist);
838 // else GetWhichTrees(f);
839
840 cout << "TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ) -- obsolte "<<endl;
841
842 SetWhichTrees(detlist);
843
844 TTree *Trout =0;
845
846 TString fname = f->GetName();
847 if(!CheckLevel2File(fname))return NULL;
848
849 // UInt_t *found=0;
850
851 cout<< "GetPamTree(TFile*,TString): detector list --> ";
852 if(TRK1)cout<<"TRK1 ";
853 if(TRK2)cout<<"TRK2 ";
854 if(TRKh)cout<<"TRKH ";
855 if(CAL1)cout<<"CAL1 ";
856 if(CAL2)cout<<"CAL2 ";
857 if(TOF)cout<<"TOF ";
858 if(TRG)cout<<"TRG ";
859 if(AC)cout<<"AC ";
860 if(ND)cout<<"ND ";
861 if(S4)cout<<"S4 ";
862 if(ORB)cout<<"ORB ";
863 cout << endl;
864 if(SELLI)cout<<">>> Found selection-list <<<"<<endl;
865
866 f->cd();
867
868 // Tracker
869 TTree *T = (TTree*)f->Get("Tracker");
870 if(T && (TRK2||TRK1||TRKh)) {
871 if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
872 // else T->SetBranchStatus("TrkLevel2",0,found);
873 if(TRK2)cout << "Tracker : set branch address TrkLevel2"<<endl;
874 if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
875 // else T->SetBranchStatus("TrkLevel1",0,found);
876 if(TRK1)cout << "Tracker : set branch address TrkLevel1"<<endl;
877 if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
878 // else T->SetBranchStatus("TrkHough",0,found);
879 if(TRKh)cout << "Tracker : set branch address TrkHough"<<endl;
880 if(!Trout)Trout=T;
881 else Trout->AddFriend(T);
882 }else{
883 cout << "Tracker : missing tree"<<endl;
884 };
885 // Calorimeter
886 TTree *C = (TTree*)f->Get("Calorimeter");
887 if(C && (CAL2||CAL1)) {
888 if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
889 // else C->SetBranchStatus("CaloLevel2",0,found);
890 if(CAL2)cout << "Calorimeter : set branch address CaloLevel2"<<endl;
891 if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
892 // else C->SetBranchStatus("CaloLevel1",0,found);
893 if(CAL1)cout << "Calorimeter : set branch address CaloLevel1"<<endl;
894 if(!Trout)Trout=C;
895 else Trout->AddFriend(C);
896 }else{
897 cout << "Calorimeter : missing tree"<<endl;
898 };
899
900 // ToF
901 TTree *O = (TTree*)f->Get("ToF");
902 if(O && TOF) {
903 O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
904 cout << "ToF : set branch address ToFLevel2"<<endl;
905 if(!Trout)Trout=O;
906 else Trout->AddFriend(O);
907 }else{
908 cout << "ToF : missing tree"<<endl;
909 };
910 // Trigger
911 TTree *R = (TTree*)f->Get("Trigger");
912 if(R && TRG) {
913 R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
914 cout << "Trigger : set branch address TrigLevel2"<<endl;
915 if(!Trout)Trout=O;
916 else Trout->AddFriend(R);
917 }else{
918 cout << "Trigger : missing tree"<<endl;
919 };
920 // S4
921 TTree *S = (TTree*)f->Get("S4");
922 if(S && S4) {
923 S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
924 cout << "S4 : set branch address S4Level2"<<endl;
925 if(!Trout)Trout=O;
926 else Trout->AddFriend(S);
927 }else{
928 cout << "S4 : missing tree"<<endl;
929 };
930 // Neutron Detector
931 TTree *N = (TTree*)f->Get("NeutronD");
932 if(N && ND) {
933 N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
934 cout << "NeutronD : set branch address NDLevel2"<<endl;
935 if(!Trout)Trout=O;
936 else Trout->AddFriend(N);
937 }else{
938 cout << "NeutronD : missing tree"<<endl;
939 };
940 // Anticounters
941 TTree *A = (TTree*)f->Get("Anticounter");
942 if(A && AC) {
943 A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
944 cout << "Anticounter : set branch address AcLevel2"<<endl;
945 if(!Trout)Trout=O;
946 else Trout->AddFriend(A);
947 }else{
948 cout << "Anticounter : missing tree"<<endl;
949 };
950 // Orbital Info
951 TTree *B = (TTree*)f->Get("OrbitalInfo");
952 if(B && ORB) {
953 B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
954 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
955 if(!Trout)Trout=O;
956 else Trout->AddFriend(B);
957 }else{
958 cout << "OrbitalInfo : missing tree"<<endl;
959 };
960
961 TTree *L = (TTree*)f->Get("SelectionList");
962 if(L && SELLI==1) {
963 // L->SetBranchAddress("RunEntry",&irun);
964 // cout << "SelectionList: set branch address RunEntry"<<endl;
965 // L->SetBranchAddress("EventEntry",&irunentry);
966 // cout << "SelectionList: set branch address EventEntry"<<endl;
967 // if(!Trout)Trout=O;
968 // else Trout->AddFriend("SelectionList");
969 cout << " TTree *PamLevel2::GetPamTree(TFile, TString) >>> SelectionList not implemente!!!!"<<endl;
970 sel_tree = 0;
971 }else{
972 cout << "SelectionList : missing tree"<<endl;
973 };
974
975 cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
976
977 pam_tree = (TChain*)Trout;
978
979 return Trout;
980
981 }
982 //--------------------------------------
983 //
984 //
985 //--------------------------------------
986 /**
987 * Get list of Level2 files.
988 * @param ddir Level2 data directory.
989 * @param flisttxt Name of txt file containing file list.
990 * @return Pointer to a TList of TSystemFiles
991 * If no input file list is given , all the Level2 files inside the directory are processed.
992 */
993 TList* PamLevel2::GetListOfLevel2Files(TString ddir, TString flisttxt = ""){
994
995 TString wdir = gSystem->WorkingDirectory();
996
997 if(ddir=="")ddir = wdir;
998 // TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
999 cout << "Level2 data directory : "<< ddir << endl;
1000
1001 TList *contents = new TList; // create output list
1002 contents->SetOwner();
1003
1004 // char *fullpath;
1005 // const char *fullpath;
1006
1007 // if no input file list is given:
1008 if ( flisttxt != "" ){
1009
1010 // if( !gSystem->IsFileInIncludePath(flisttxt,&fullpath) ){
1011 // if( !(fullpath = gSystem->FindFile("./",flisttxt)) ){
1012 // cout <<"File "<<flisttxt<<" not found"<<endl;
1013 // return 0;
1014 // }
1015 // flisttxt = fullpath;
1016 if ( !flisttxt.EndsWith(".root") ){
1017
1018 flisttxt = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1019
1020 if( !gSystem->ChangeDirectory(ddir) ){
1021 cout << "Cannot change directory : "<<ddir<<endl;
1022 return 0;
1023 }
1024
1025 cout <<"Input file list : " << flisttxt <<endl;
1026 ifstream in;
1027 in.open(flisttxt, ios::in); //open input file list
1028 if(!in.good()){
1029 cout <<" ERROR opening the file "<<endl;
1030 gSystem->ChangeDirectory(wdir); // back to the working directory
1031 return 0;
1032 }
1033 int line=0;
1034 while (1) {
1035 TString file;
1036 in >> file;
1037 if (!in.good()) break;
1038 line++;
1039 // cout <<"(1) " << file << endl;
1040 if(file.IsNull()){
1041 cout << "-- list interrupted at line "<<line <<endl;
1042 break;
1043 }
1044 if(file.Contains("#"))file = file(0,file.First("#"));
1045 // cout <<"(2) " << file << endl;
1046 // if( gSystem->IsFileInIncludePath(file,&fullpath) ){
1047 // if( (fullpath = gSystem->FindFile(ddir,file)) ){
1048 if( file.EndsWith(".root") ){
1049 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(ddir),gSystem->BaseName(file));
1050 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1051 delete fullpath;
1052 }
1053 // }else{
1054 // if(file.Data()!="")cout << "File: "<<file<<" ---> missing "<< endl;
1055 // };
1056 };
1057 in.close();
1058 } else {
1059 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1060 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1061 delete fullpath;
1062 };
1063 }else{
1064
1065 cout << "No input file list given."<<endl;
1066 cout << "Check for existing root files."<<endl;
1067 // cout << "Warking directory: "<< gSystem->WorkingDirectory()<< endl;
1068
1069 TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
1070 TList *temp = datadir->GetListOfFiles();
1071 // temp->Print();
1072 // cout << "*************************************" << endl;
1073
1074 TIter next(temp);
1075 TSystemFile *questo = 0;
1076
1077
1078 while ( (questo = (TSystemFile*) next()) ) {
1079 TString name = questo-> GetName();
1080 if( name.EndsWith(".root") ){
1081 // const char *fullpath = gSystem->FindFile(ddir,name);
1082 // char *fullpath;
1083 // gSystem->IsFileInIncludePath(name,&fullpath);
1084 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(ddir),gSystem->BaseName(name));
1085 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));
1086 delete fullpath;
1087 };
1088 }
1089 delete temp;
1090 delete datadir;
1091
1092 };
1093 gSystem->ChangeDirectory(wdir); // back to the working directory
1094 // cout << endl << "Selected files:" << endl;
1095 // contents->Print();
1096 cout << contents->GetEntries()<<" files \n";
1097 // cout << endl;
1098 // cout << "Working directory: "<< gSystem->WorkingDirectory()<< endl;
1099 return contents;
1100 };
1101 //--------------------------------------
1102 //
1103 //
1104 //--------------------------------------
1105 /**
1106 * Get the Pamela detector chains from a list of files and make them friends.
1107 * @param fl Pointer to a TList of TSystemFiles
1108 * @param detlist String to select trees to be included
1109 * @return Pointer to a TChain
1110 */
1111 TChain *PamLevel2::GetPamTree(TList *fl, TString detlist ){
1112
1113 TChain *Trout =0;
1114
1115 // if( !detlist.IsNull() )SetWhichTrees(detlist);
1116 SetWhichTrees(detlist);
1117
1118 cout<< "GetPamTree(TList*,TString): input detector list --> ";
1119 if(TRK1)cout<<"TRK1 ";
1120 if(TRK2)cout<<"TRK2 ";
1121 if(TRKh)cout<<"TRKH ";
1122 if(CAL1)cout<<"CAL1 ";
1123 if(CAL2)cout<<"CAL2 ";
1124 if(TOF)cout<<"TOF ";
1125 if(TRG)cout<<"TRG ";
1126 if(AC)cout<<"AC ";
1127 if(ND)cout<<"ND ";
1128 if(S4)cout<<"S4 ";
1129 if(ORB)cout<<"ORB ";
1130 cout << endl;
1131
1132 TChain *T = 0;
1133 TChain *C = 0;
1134 TChain *O = 0;
1135 TChain *R = 0;
1136 TChain *S = 0;
1137 TChain *N = 0;
1138 TChain *A = 0;
1139 TChain *B = 0;
1140
1141 TChain *L = 0;
1142
1143 if(TRK2||TRK1||TRKh) T = new TChain("Tracker");
1144 if(CAL2||CAL1) C = new TChain("Calorimeter");
1145 if(TOF) O = new TChain("ToF");
1146 if(TRG) R = new TChain("Trigger");
1147 if(S4) S = new TChain("S4");
1148 if(ND) N = new TChain("NeutronD");
1149 if(AC) A = new TChain("Anticounter");
1150 if(ORB) B = new TChain("OrbitalInfo");
1151
1152 L = new TChain("SelectionList");
1153
1154 // loop over files and create chains
1155 TIter next(fl);
1156 TSystemFile *questo = 0;
1157 while ( (questo = (TSystemFile*) next()) ) {
1158 TString name = questo->GetName();
1159 cout << "File: "<< name << endl;
1160 if( CheckLevel2File(name) ){
1161 if(TRK2||TRK1||TRKh) T->Add(name);
1162 if(CAL1||CAL2) C->Add(name);
1163 if(TOF) O->Add(name);
1164 if(TRG) R->Add(name);
1165 if(S4) S->Add(name);
1166 if(ND) N->Add(name);
1167 if(AC) A->Add(name);
1168 if(ORB) B->Add(name);
1169 if(SELLI==1) L->Add(name);
1170 };
1171 }
1172
1173 cout << "done chain \n";
1174
1175 // UInt_t *found=0;
1176 // Tracker
1177 if(T && (TRK2||TRK1||TRKh)) {
1178 if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1179 // else T->SetBranchStatus("TrkLevel2",0,found);
1180 if(TRK2)cout << "Tracker : set branch address TrkLevel2"<<endl;
1181 if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1182 // else T->SetBranchStatus("TrkLevel1",0,found);
1183 if(TRK1)cout << "Tracker : set branch address TrkLevel1"<<endl;
1184 if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1185 // else T->SetBranchStatus("TrkHough",0,found);
1186 if(TRKh)cout << "Tracker : set branch address TrkHough"<<endl;
1187 if(!Trout)Trout=T;
1188 else Trout->AddFriend("Tracker");
1189 }else{
1190 cout << "Tracker : missing tree"<<endl;
1191 };
1192 // Calorimeter
1193 if(C && (CAL2||CAL1)) {
1194 if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1195 // else C->SetBranchStatus("CaloLevel2",0,found);
1196 if(CAL2)cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1197 if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1198 // else C->SetBranchStatus("CaloLevel1",0,found);
1199 if(CAL1)cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1200 if(!Trout)Trout=C;
1201 else Trout->AddFriend("Calorimeter");
1202 }else{
1203 cout << "Calorimeter : missing tree"<<endl;
1204 };
1205 // ToF
1206 if(O && TOF) {
1207 O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1208 cout << "ToF : set branch address ToFLevel2"<<endl;
1209 if(!Trout)Trout=O;
1210 else Trout->AddFriend("ToF");
1211 }else{
1212 cout << "ToF : missing tree"<<endl;
1213 };
1214 // Trigger
1215 if(R && TRG) {
1216 R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1217 cout << "Trigger : set branch address TrigLevel2"<<endl;
1218 if(!Trout)Trout=O;
1219 else Trout->AddFriend("Trigger");
1220 }else{
1221 cout << "Trigger : missing tree"<<endl;
1222 };
1223 // S4
1224 if(S && S4) {
1225 S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1226 cout << "S4 : set branch address S4Level2"<<endl;
1227 if(!Trout)Trout=O;
1228 else Trout->AddFriend("S4");
1229 }else{
1230 cout << "S4 : missing tree"<<endl;
1231 };
1232 // Neutron Detector
1233 if(N && ND) {
1234 N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1235 cout << "NeutronD : set branch address NDLevel2"<<endl;
1236 if(!Trout)Trout=O;
1237 else Trout->AddFriend("NeutronD");
1238 }else{
1239 cout << "NeutronD : missing tree"<<endl;
1240 };
1241 // Anticounters
1242 if(A && AC) {
1243 A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1244 cout << "Anticounter : set branch address AcLevel2"<<endl;
1245 if(!Trout)Trout=O;
1246 else Trout->AddFriend("Anticounter");
1247 }else{
1248 cout << "Anticounter : missing tree"<<endl;
1249 };
1250 // Orbital Info
1251 if(B && ORB) {
1252 B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1253 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1254 if(!Trout)Trout=O;
1255 else Trout->AddFriend("OrbitalInfo");
1256 }else{
1257 cout << "OrbitalInfo : missing tree"<<endl;
1258 };
1259
1260 // Selection List
1261 if(L && SELLI==1) {
1262 cout<<">>> Found selection-list <<<"<<endl;
1263 L->SetBranchAddress("RunEntry",&irun);
1264 cout << "SelectionList: set branch address RunEntry"<<endl;
1265 L->SetBranchAddress("EventEntry",&irunentry);
1266 cout << "SelectionList: set branch address EventEntry"<<endl;
1267 sel_tree = L;
1268 // if(!Trout)Trout=O;
1269 // else Trout->AddFriend("SelectionList");
1270 }else{
1271 // cout << "SelectionList : missing tree"<<endl;
1272 if(L)L->Delete();
1273 };
1274
1275 // cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
1276
1277 pam_tree = Trout;
1278
1279 return Trout;
1280 }
1281
1282
1283
1284 //--------------------------------------
1285 //
1286 //
1287 //--------------------------------------
1288 /**
1289 * Set branch addresses for Pamela friend trees
1290 */
1291 void PamLevel2::SetBranchAddress(TTree *t){
1292
1293 TRK2 = TRK2 & t->GetBranchStatus("TrkLevel2");
1294 TRK1 = TRK1 & t->GetBranchStatus("TrkLevel1");
1295 TRKh = TRKh & t->GetBranchStatus("TrkHough");
1296 CAL2 = CAL2 & t->GetBranchStatus("CaloLevel2");
1297 CAL1 = CAL1 & t->GetBranchStatus("CaloLevel1");
1298 TOF = TOF & t->GetBranchStatus("ToFLevel2");
1299 TRG = TRG & t->GetBranchStatus("TrigLevel2");
1300 S4 = S4 & t->GetBranchStatus("S4Level2");
1301 ND = ND & t->GetBranchStatus("NDLevel2");
1302 AC = AC & t->GetBranchStatus("AcLevel2");
1303 ORB = ORB & t->GetBranchStatus("OrbitalInfo");
1304
1305
1306 // Tracker
1307 if(TRK1) {
1308 t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel2"));
1309 cout << "Tracker : set branch address TrkLevel1"<<endl;
1310 };
1311 if(TRK2) {
1312 t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel1"));
1313 cout << "Tracker : set branch address TrkLevel2"<<endl;
1314 };
1315 if(TRKh) {
1316 t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1317 cout << "Tracker : set branch address TrkHough"<<endl;
1318 };
1319
1320 // Calorimeter
1321 if(CAL1) {
1322 t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1323 cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1324 };
1325 if(CAL2) {
1326 t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1327 cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1328 };
1329
1330 // ToF
1331 if(TOF) {
1332 t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1333 cout << "ToF : set branch address ToFLevel2"<<endl;
1334 };
1335 // Trigger
1336 if(TRG) {
1337 t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1338 cout << "Trigger : set branch address TrigLevel2"<<endl;
1339 };
1340 // S4
1341 if(S4) {
1342 t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1343 cout << "S4 : set branch address S4Level2"<<endl;
1344 };
1345 // Neutron Detector
1346 if(ND) {
1347 t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1348 cout << "NeutronD : set branch address NDLevel2"<<endl;
1349 };
1350 // Anticounters
1351 if(AC) {
1352 t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1353 cout << "Anticounter : set branch address AcLevel2"<<endl;
1354 };
1355 // OrbitalInfo
1356 if(ORB) {
1357 t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1358 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1359 };
1360 // SelectionList
1361 if(SELLI==1) {
1362 t->SetBranchAddress("RunEntry", &irun);
1363 cout << "SelectionList: set branch address RunEntry"<<endl;
1364 t->SetBranchAddress("EventEntry", &irunentry);
1365 cout << "SelectionList: set branch address EventEntry"<<endl;
1366 };
1367
1368 }
1369 /**
1370 * Set branch addresses for Pamela friend trees
1371 */
1372 void PamLevel2::SetBranchAddress(TChain *t){
1373
1374 TRK2 = TRK2 & t->GetBranchStatus("TrkLevel2");
1375 TRK1 = TRK1 & t->GetBranchStatus("TrkLevel1");
1376 TRKh = TRKh & t->GetBranchStatus("TrkHough");
1377 CAL1 = CAL1 & t->GetBranchStatus("CaloLevel1");
1378 CAL2 = CAL2 & t->GetBranchStatus("CaloLevel2");
1379 TOF = TOF & t->GetBranchStatus("ToFLevel2");
1380 TRG = TRG & t->GetBranchStatus("TrigLevel2");
1381 S4 = S4 & t->GetBranchStatus("S4Level2");
1382 ND = ND & t->GetBranchStatus("NDLevel2");
1383 AC = AC & t->GetBranchStatus("AcLevel2");
1384 ORB = ORB & t->GetBranchStatus("OrbitalInfo");
1385
1386 // Tracker
1387 if(TRK2) {
1388 t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1389 cout << "Tracker : set branch address TrkLevel2"<<endl;
1390 };
1391 if(TRK1) {
1392 t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1393 cout << "Tracker : set branch address TrkLevel1"<<endl;
1394 };
1395 if(TRKh) {
1396 t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1397 cout << "Tracker : set branch address TrkHough"<<endl;
1398 };
1399
1400 // Calorimeter
1401 if(CAL2) {
1402 t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1403 cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1404 };
1405 if(CAL1) {
1406 t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1407 cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1408 };
1409
1410 // ToF
1411 if(TOF) {
1412 t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1413 cout << "ToF : set branch address ToFLevel2"<<endl;
1414 };
1415 // Trigger
1416 if(TRG) {
1417 t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1418 cout << "Trigger : set branch address TrigLevel2"<<endl;
1419 };
1420 // S4
1421 if(S4) {
1422 t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1423 cout << "S4 : set branch address S4Level2"<<endl;
1424 };
1425 // Neutron Detector
1426 if(ND) {
1427 t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1428 cout << "NeutronD : set branch address NDLevel2"<<endl;
1429 };
1430 // Anticounters
1431 if(AC) {
1432 t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1433 cout << "Anticounter : set branch address AcLevel2"<<endl;
1434 };
1435 // OrbitalInfo
1436 if(ORB) {
1437 t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1438 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1439 };
1440 // SelectionList
1441 if(SELLI==1) {
1442 t->SetBranchAddress("RunEntry", &irun);
1443 cout << "SelectionList: set branch address RunEntry"<<endl;
1444 t->SetBranchAddress("EventEntry", &irunentry);
1445 cout << "SelectionList: set branch address EventEntry"<<endl;
1446 };
1447
1448 }
1449
1450
1451 //--------------------------------------
1452 //
1453 //
1454 //--------------------------------------
1455 /**
1456 * Get the Run tree chain from a list of files.
1457 * @param fl Pointer to a TList of TSystemFiles
1458 * @return Pointer to a TChain
1459 */
1460 TChain *PamLevel2::GetRunTree(TList *fl){
1461
1462 TChain *R = new TChain("Run");
1463
1464 // loop over files and create chains
1465 TIter next(fl);
1466 TSystemFile *questo = 0;
1467 while ( (questo = (TSystemFile*) next()) ) {
1468 TString name = questo->GetName();
1469 // cout << "File: "<< name << endl;
1470 if( CheckLevel2File(name) ){
1471 R->Add(name);
1472 };
1473 }
1474
1475 if(R->GetNtrees()){
1476
1477 R->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1478 cout << "Run : set branch address RunInfo"<<endl;
1479 R->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1480 cout << "Software : set branch address SoftInfo"<<endl; // Emiliano
1481 }else{
1482 delete R;
1483 R=0;
1484 }
1485
1486 run_tree = R;
1487
1488 return R;
1489
1490 }
1491 //--------------------------------------
1492 //
1493 //
1494 //--------------------------------------
1495 /**
1496 * Get the Run tree from a file.
1497 * @param f Pointer to a TFile
1498 * @return Pointer to a TTree
1499 */
1500 TTree *PamLevel2::GetRunTree(TFile *f){
1501
1502 cout << "TTree *PamLevel2::GetRunTree(TFile *f) -- obsolte "<<endl;
1503
1504 TTree *R = (TTree*)f->Get("Run");
1505
1506 if(R){
1507 R->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1508 cout << "Run : set branch address RunInfo"<<endl;
1509 R->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1510 cout << "Software : set branch address SoftInfo"<<endl; // Emiliano
1511 }
1512
1513 run_tree = (TChain*)R;
1514
1515 return R;
1516
1517 }
1518 /**
1519 * Update the runinfo informations (to be used to have Run infos event by event basis)
1520 * @param run Pointer to the chain/tree which contains run infos
1521 * @return true if a new run has been read, false if it is still the same run
1522 */
1523 Bool_t PamLevel2::UpdateRunInfo(TChain *run, Long64_t iev){
1524 //
1525 // check if we have already called once GetEntry, if not call it
1526 //
1527 if(!run){
1528 cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing RunInfo tree "<<endl;
1529 return(false);
1530 }
1531 if ( run->GetEntries() <= 0 ) return(false);
1532 //
1533
1534 // Int_t oldrun = irun;
1535 Long64_t oldrun = irun;
1536 // --------------------------------------
1537 // if it is a full file (not preselected)
1538 // --------------------------------------
1539 if(SELLI==0){
1540
1541 if( !GetOrbitalInfo() ){
1542 cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
1543 return(false);
1544 }
1545
1546 if ( irun < 0 ){
1547 irun = 0;
1548 run->GetEntry(irun);
1549 runfirstentry = 0LL;
1550 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1551 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1552 // cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1553 // cout << "runfirstentry "<<runfirstentry<<endl;
1554 // return(true);
1555 };
1556
1557 // modificato il controllo sull'aggiornamento del run, per evitare problemi
1558 // dovuti agli eventi annidati (NB! NEVENTS conta anche questi!!)
1559 //
1560 bool fromfirst = true;
1561 //
1562 while ( !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) && irun < run->GetEntries()-1 ){
1563 // while ( iev > (runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS-1)) && irun < run->GetEntries() ){
1564 irun++;
1565 run->GetEntry(irun);
1566 runfirstentry = runlastentry;
1567 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runfirstentry += 1LL;
1568 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1569 // cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1570 // cout << "runfirstentry "<<runfirstentry<<endl;
1571 // printf(" iev %llu %u %llu \n",iev,this->GetRunInfo()->NEVENTS,(ULong64_t)(runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS)));
1572 // printf(" abstime %u trailertime %u \n",GetOrbitalInfo()->absTime,GetRunInfo()->RUNTRAILER_TIME);
1573 // printf(" IDRUN %u \n",GetRunInfo()->ID);
1574 //
1575 // prevshift = 0;
1576 //
1577 if ( irun == run->GetEntries()-1 && fromfirst && !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)){
1578 printf(" resetting irun (it should NOT happen!!!)\n");
1579 fromfirst = false;
1580 irun = 0;
1581 run->GetEntry(irun);
1582 runfirstentry = 0ULL;
1583 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1584 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1585 };
1586 //
1587 };
1588 //
1589 if ( !fromfirst &&
1590 !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&
1591 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)
1592 ) {
1593 printf(" Something very wrong here: cannot find RUN containing absolute time %u \n",GetOrbitalInfo()->absTime);
1594 return false;
1595 }
1596 //
1597 if ( irun == oldrun || irun >= run->GetEntries() ) return(false);
1598 //
1599 // printf(" iev %llu irun %i nevents %u 1st %llu last %llu \n",iev,irun,this->GetRunInfo()->NEVENTS,(ULong64_t)runfirstentry,(ULong64_t)runlastentry);
1600 //
1601 prevshift = 0;
1602 cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1603 // cout << "runfirstentry "<<runfirstentry<<endl;
1604 return(true);
1605 };
1606 // ----------------------------------------------------
1607 // if it is a preselected file (there is SelectionList)
1608 // NBNB - the event tree MUST be read first
1609 // ----------------------------------------------------
1610 if(SELLI==1){
1611 sel_tree->GetEntry(iev);
1612 // cout << irun << " "<< irunentry << endl;
1613 if(irun != oldrun){
1614 run->GetEntry(irun);
1615 cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1616 prevshift = 0;
1617 return true;
1618 }
1619 return false;
1620 }
1621
1622 return false;
1623 //
1624 };
1625
1626 // Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev){
1627 // //
1628 // // check if we have already called once GetEntry, if not call it
1629 // //
1630 // if(!run_tree){
1631 // cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree not loaded"<<endl;
1632 // return false;
1633 // }
1634 // if ( run_tree->GetEntries() <= 0 ) {
1635 // cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree is empty"<<endl;
1636 // return(false);
1637 // }
1638 // //
1639
1640 // Int_t oldrun = irun;
1641 // // --------------------------------------
1642 // // if it is a full file (not preselected)
1643 // // --------------------------------------
1644 // if(SELLI==0){
1645
1646 // if( !GetOrbitalInfo() ){
1647 // cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
1648 // return(false);
1649 // }
1650
1651 // if ( irun < 0 ){
1652 // irun = 0;
1653 // run_tree->GetEntry(irun);
1654 // runfirstentry = iev;
1655 // if( iev != 0){
1656 // cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- WARNING -- "<<endl;
1657 // }
1658 // };
1659
1660 // bool fromfirst = true;
1661
1662 // while ( !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) && irun < run->GetEntries() ){
1663 // irun++;
1664 // if( irun == run->GetEntries() ){
1665 // if(!fromfirst){
1666 // cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- event entry #"<<iev<<" does not belong to any run (should not happen)" <<endl;
1667 // return false;
1668 // }
1669 // cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- reached end of run tree. searchin again from beginning (should not happen)" <<endl;
1670 // fromfirst = false;
1671 // irun =0;
1672 // }
1673 // run_tree->GetEntry(irun);
1674 // runfirstentry = iev;
1675 // };
1676 // //
1677 // if ( irun == oldrun ) return(false);
1678 // //
1679 // cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1680 // cout << "runfirstentry "<<runfirstentry<<endl;
1681 // prevshift = 0;
1682 // return(true);
1683 // };
1684 // // ----------------------------------------------------
1685 // // if it is a preselected file (there is SelectionList)
1686 // // NBNB - the event tree MUST be read first
1687 // // ----------------------------------------------------
1688 // if(SELLI==1){
1689 // sel_tree->GetEntry(iev);
1690 // if(irun != oldrun){
1691 // run->GetEntry(irun);
1692 // cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1693 // prevshift = 0;
1694 // return true;
1695 // }
1696 // return false;
1697 // }
1698
1699 // return false;
1700 // //
1701 // };
1702 /**
1703 * Update the runinfo informations (to be used to have Run infos event by event basis)
1704 * @param run Pointer to the chain/tree which contains run infos
1705 * @return true if a new run has been read, false if it is still the same run
1706 */
1707 Bool_t PamLevel2::UpdateRunInfo(TTree *run, Long64_t iev){
1708 return(UpdateRunInfo((TChain*)run,iev));
1709 };
1710
1711 //--------------------------------------
1712 //
1713 //
1714 //--------------------------------------
1715 /**
1716 * Set which trees shoul be analysed
1717 * @param detlist TString containing the sequence of trees required
1718 */
1719 void PamLevel2::SetWhichTrees(TString detlist){
1720
1721 if(detlist.IsNull() || detlist.Contains("+ALL", TString::kIgnoreCase)){
1722 CAL0 = false;
1723 CAL1 = true;
1724 CAL2 = true;
1725 TRK2 = true;
1726 TRK1 = false;
1727 TRKh = false;
1728 TRK0 = false;
1729 TRG = true;
1730 TOF = true;
1731 TOF0 = false;
1732 S4 = true;
1733 ND = true;
1734 AC = true;
1735 ORB = true;
1736 }else if( detlist.Contains("-ALL", TString::kIgnoreCase) ){
1737 CAL0 = false;
1738 CAL1 = false;
1739 CAL2 = false;
1740 TRK2 = false;
1741 TRK1 = false;
1742 TRKh = false;
1743 TRK0 = false;
1744 TRG = false;
1745 TOF = false;
1746 TOF0 = false;
1747 S4 = false;
1748 ND = false;
1749 AC = false;
1750 ORB = false;
1751 };
1752
1753 // -------------------------------------------------------------------------
1754 if( detlist.Contains("CAL1", TString::kIgnoreCase) ){
1755 if ( detlist.Contains("-CAL1", TString::kIgnoreCase) )CAL1=false;
1756 if ( detlist.Contains("+CAL1", TString::kIgnoreCase) )CAL1=true;
1757 };
1758
1759 if( detlist.Contains("CAL0", TString::kIgnoreCase) ){
1760 if ( detlist.Contains("-CAL0", TString::kIgnoreCase) )CAL0=false;
1761 if ( detlist.Contains("+CAL0", TString::kIgnoreCase) )CAL0=true;
1762 };
1763
1764 if( detlist.Contains("CAL2", TString::kIgnoreCase)){
1765 if ( detlist.Contains("-CAL2", TString::kIgnoreCase) )CAL2=false;
1766 if ( detlist.Contains("+CAL2", TString::kIgnoreCase) )CAL2=true;
1767 };
1768
1769 if( detlist.Contains("+CAL", TString::kIgnoreCase) && !CAL1 && !CAL2 )CAL2=true;
1770 if( detlist.Contains("-CAL", TString::kIgnoreCase) && CAL1 && CAL2 ){
1771 CAL2=false;
1772 CAL1=false;
1773 }
1774 // -------------------------------------------------------------------------
1775 if( detlist.Contains("TRK0", TString::kIgnoreCase) ){
1776 if ( detlist.Contains("-TRK0", TString::kIgnoreCase) )TRK0=false;
1777 if ( detlist.Contains("+TRK0", TString::kIgnoreCase) )TRK0=true;
1778 };
1779
1780 if( detlist.Contains("TRK1", TString::kIgnoreCase) ){
1781 if ( detlist.Contains("-TRK1", TString::kIgnoreCase) )TRK1=false;
1782 if ( detlist.Contains("+TRK1", TString::kIgnoreCase) )TRK1=true;
1783 };
1784
1785 if( detlist.Contains("TRK2", TString::kIgnoreCase)){
1786 if ( detlist.Contains("-TRK2", TString::kIgnoreCase) )TRK2=false;
1787 if ( detlist.Contains("+TRK2", TString::kIgnoreCase) )TRK2=true;
1788 };
1789
1790 if( detlist.Contains("TRKh", TString::kIgnoreCase)){
1791 if ( detlist.Contains("-TRKh", TString::kIgnoreCase) )TRKh=false;
1792 if ( detlist.Contains("+TRKh", TString::kIgnoreCase) )TRKh=true;
1793 };
1794
1795 if( detlist.Contains("+TRK", TString::kIgnoreCase) && !TRK1 && !TRK2 && !TRKh )TRK2=true;
1796 if( detlist.Contains("-TRK", TString::kIgnoreCase) && TRK1 && TRK2 && TRKh){
1797 TRK2=false;
1798 TRK1=false;
1799 TRKh=false;
1800 }
1801 // -------------------------------------------------------------------------
1802
1803 if( detlist.Contains("-TRG", TString::kIgnoreCase) )TRG = false;
1804 else if( detlist.Contains("+TRG", TString::kIgnoreCase) )TRG = true;
1805
1806 if( detlist.Contains("-TOF", TString::kIgnoreCase) )TOF = false;
1807 else if( detlist.Contains("+TOF", TString::kIgnoreCase) )TOF = true;
1808
1809 if( detlist.Contains("-TOF0", TString::kIgnoreCase) )TOF0 = false;
1810 else if( detlist.Contains("+TOF0", TString::kIgnoreCase) )TOF0 = true;
1811
1812 if( detlist.Contains("-S4", TString::kIgnoreCase) )S4 = false;
1813 else if( detlist.Contains("+S4", TString::kIgnoreCase) )S4 = true;
1814
1815 if( detlist.Contains("-ND", TString::kIgnoreCase) )ND = false;
1816 else if( detlist.Contains("+ND", TString::kIgnoreCase) )ND = true;
1817
1818 if( detlist.Contains("-AC", TString::kIgnoreCase) )AC = false;
1819 else if( detlist.Contains("+AC", TString::kIgnoreCase) )AC = true;
1820
1821 if( detlist.Contains("-ORB", TString::kIgnoreCase) )ORB = false;
1822 else if( detlist.Contains("+ORB", TString::kIgnoreCase) )ORB = true;
1823
1824 // cout<< "Set detector list --> ";
1825 // if(TRK1)cout<<"TRK1 ";
1826 // if(TRK2)cout<<"TRK2 ";
1827 // if(TRKh)cout<<"TRKH ";
1828 // if(CAL1)cout<<"CAL1 ";
1829 // if(CAL2)cout<<"CAL2 ";
1830 // if(TOF)cout<<"TOF ";
1831 // if(TRG)cout<<"TRG ";
1832 // if(AC)cout<<"AC ";
1833 // if(ND)cout<<"ND ";
1834 // if(S4)cout<<"S4 ";
1835 // if(ORB)cout<<"ORB ";
1836 // cout << endl;
1837
1838 };
1839
1840
1841 /**
1842 * Set tree/branch detector flags from the content of a tree
1843 */
1844 void PamLevel2::GetWhichTrees(TFile* f){
1845
1846
1847
1848 cout << "void PamLevel2::GetWhichTrees(TFile* f) --- WARNING!! --- ...potrebbe non funzionare "<<endl;
1849 // -----------
1850 // reset flags
1851 // -----------
1852 CAL1 = false;
1853 CAL2 = false;
1854 TRK2 = false;
1855 TRK1 = false;
1856 TRKh = false;
1857 TRG = false;
1858 TOF = false;
1859 S4 = false;
1860 ND = false;
1861 AC = false;
1862 ORB = false;
1863
1864 RUN = false;
1865
1866 cout << "Checking file: "<<f->GetName()<<endl;
1867 if( !f || f->IsZombie() ){
1868 cout << "File: "<< f->GetName() <<" Non valid root file"<< endl;
1869 return;
1870 }
1871
1872 TList *lk = f->GetListOfKeys();
1873 TIter next(lk);
1874 TKey *key =0;
1875
1876 Int_t nev = 0;
1877
1878 while( (key = (TKey*)next()) ){
1879
1880 if( !strcmp(key->GetName(),"Run" ) )RUN = true;
1881
1882 //=========================================================
1883 if( !strcmp(key->GetName(),"Trigger" ) ){
1884 TRG = true;
1885 Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
1886 if( nev && nevt!=nev){
1887 cout << "File: "<< f->GetName() <<" Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
1888 TRG = false;
1889 }else nev=nevt;
1890 }
1891 //=========================================================
1892 if( !strcmp(key->GetName(),"ToF" ) ){
1893 TOF = true;
1894 Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
1895 if( nev && nevt!=nev){
1896 cout << "File: "<< f->GetName() <<" ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
1897 TOF = false;
1898 }else nev=nevt;
1899 }
1900 //=========================================================
1901 if( !strcmp(key->GetName(),"S4" ) ){
1902 S4 = true;
1903 Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
1904 if( nev && nevt!=nev){
1905 cout << "File: "<< f->GetName() <<" S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
1906 S4 = false;
1907 }else nev=nevt;
1908 }
1909 //=========================================================
1910
1911 if( !strcmp(key->GetName(),"NeutronD" ) ){
1912 ND = true;
1913 Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
1914 if( nev && nevt!=nev){
1915 cout << "File: "<< f->GetName() <<"NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
1916 ND =false;
1917 }else nev=nevt;
1918 }
1919 //=========================================================
1920 if( !strcmp(key->GetName(),"Anticounter") ){
1921 AC = true;
1922 Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
1923 if( nev && nevt!=nev){
1924 cout << "File: "<< f->GetName() <<" Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
1925 AC =false;
1926 }else nev=nevt;
1927 }
1928 //=========================================================
1929 if( !strcmp(key->GetName(),"OrbitalInfo") ){
1930 ORB = true;
1931 Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
1932 if( nev && nevt!=nev){
1933 cout << "File: "<< f->GetName() <<" OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
1934 ORB = false;
1935 }else nev=nevt;
1936 }
1937 //=========================================================
1938 if( !strcmp(key->GetName(),"Tracker" ) ){
1939 TTree *T = (TTree*)f->Get("Tracker");
1940 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
1941 TString name = T->GetListOfBranches()->At(i)->GetName();
1942 if( !name.CompareTo("TrkLevel1") )TRK1=true;
1943 if( !name.CompareTo("TrkLevel2") )TRK2=true;
1944 if( !name.CompareTo("TrkHough") )TRKh=true;
1945 };
1946 Int_t nevt = T->GetEntries();
1947 if( nev && nevt!=nev){
1948 cout << "File: "<< f->GetName() <<" Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
1949 TRK1 = false;
1950 TRK2 = false;
1951 TRKh = false;
1952 }else nev=nevt;
1953 T->Delete();
1954 };
1955 //=========================================================
1956 if( !strcmp(key->GetName(),"Calorimeter" ) ){
1957 TTree *T = (TTree*)f->Get("Calorimeter");
1958 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
1959 TString name = T->GetListOfBranches()->At(i)->GetName();
1960 if( !name.CompareTo("CaloLevel1") )CAL1=true;
1961 if( !name.CompareTo("CaloLevel2") )CAL2=true;
1962 };
1963 Int_t nevt = T->GetEntries();
1964 if( nev && nevt!=nev){
1965 cout << "File: "<< f->GetName() <<" Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
1966 CAL1 = false;
1967 CAL2 = false;
1968 }else nev=nevt;
1969 T->Delete();
1970 };
1971
1972 };
1973
1974 delete lk;
1975
1976 // cout<< "Get detector list from input file --> ";
1977 // if(TRK1)cout<<"TRK1 ";
1978 // if(TRK2)cout<<"TRK2 ";
1979 // if(TRKh)cout<<"TRKH ";
1980 // if(CAL1)cout<<"CAL1 ";
1981 // if(CAL2)cout<<"CAL2 ";
1982 // if(TOF)cout<<"TOF ";
1983 // if(TRG)cout<<"TRG ";
1984 // if(AC)cout<<"AC ";
1985 // if(ND)cout<<"ND ";
1986 // if(S4)cout<<"S4 ";
1987 // if(ORB)cout<<"ORB ";
1988 // cout << endl;
1989
1990 return ;
1991
1992 };
1993
1994
1995 //--------------------------------------
1996 //
1997 //
1998 //--------------------------------------
1999 /**
2000 * Check if a file contains selected Pamela Level2 trees.
2001 * @param name File name
2002 * @return true if the file is ok.
2003 */
2004 Bool_t PamLevel2::CheckLevel2File(TString name){
2005
2006 Bool_t CAL1__ok = false;
2007 Bool_t CAL2__ok = false;
2008 Bool_t TRK2__ok = false;
2009 Bool_t TRK1__ok = false;
2010 Bool_t TRKh__ok = false;
2011 Bool_t TRG__ok = false;
2012 Bool_t TOF__ok = false;
2013 Bool_t S4__ok = false;
2014 Bool_t ND__ok = false;
2015 Bool_t AC__ok = false;
2016 Bool_t ORB__ok = false;
2017
2018 Bool_t RUN__ok = false;
2019
2020 Bool_t SELLI__ok = false;
2021
2022 cout << "Checking file: "<<name<<endl;
2023 TFile *f = new TFile(name.Data());
2024 if( !f || f->IsZombie() ){
2025 cout << "File: "<< f->GetName() <<" discarded ---- Non valid root file"<< endl; return false;
2026 }
2027 // cout << "Get list of keys: "<<f<<endl;
2028 TList *lk = f->GetListOfKeys();
2029 // lk->Print();
2030 TIter next(lk);
2031 TKey *key =0;
2032
2033 Int_t nev = 0;
2034
2035 while( (key = (TKey*)next()) ){
2036
2037 // cout << key->GetName() << endl;
2038 // cout << key->GetName() << ""<<key->GetClassName()<<endl;
2039 // cout << " Get tree: " << f->Get(key->GetName())<<endl;
2040 // nev_previous = nev;
2041 // cout << " n.entries "<< nev <<endl;
2042 // if( key->GetClassName()=="TTree" && nev_previous && nev != nev_previous ){
2043 // nev = ((TTree*)f->Get(key->GetName()))->GetEntries();
2044 // cout << "File: "<< f->GetName() <<" discarded ---- "<< key->GetName() << " tree: n.entries does not match "<<nev<<" "<<nev_previous<< endl;
2045 // return false;
2046 // };
2047
2048
2049 if( !strcmp(key->GetName(),"Run" ) )RUN__ok = true;
2050
2051 //=========================================================
2052 if( !strcmp(key->GetName(),"SelectionList" ) ){
2053 SELLI__ok = true;
2054 if(SELLI==1){
2055 Int_t nevt = ((TTree*)f->Get("SelectionList"))->GetEntries();
2056 if( nev && nevt!=nev){
2057 cout << "File: "<< f->GetName() <<" discarded ---- SelectionList tree has "<<nevt<<" events instead of "<<nev<< endl;
2058 return false;
2059 }
2060 nev=nevt;
2061 }
2062 }
2063
2064 //=========================================================
2065 if( !strcmp(key->GetName(),"Trigger" ) ){
2066 TRG__ok = true;
2067 if(TRG){
2068 Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
2069 if( nev && nevt!=nev){
2070 cout << "File: "<< f->GetName() <<" discarded ---- Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
2071 return false;
2072 }
2073 nev=nevt;
2074 }
2075 }
2076 //=========================================================
2077 if( !strcmp(key->GetName(),"ToF" ) ){
2078 TOF__ok = true;
2079 if(TOF){
2080 Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
2081 if( nev && nevt!=nev){
2082 cout << "File: "<< f->GetName() <<" discarded ---- ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
2083 return false;
2084 }
2085 nev=nevt;
2086 }
2087 }
2088 //=========================================================
2089 if( !strcmp(key->GetName(),"S4" ) ){
2090 S4__ok = true;
2091 if(S4){
2092 Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
2093 if( nev && nevt!=nev){
2094 cout << "File: "<< f->GetName() <<" discarded ---- S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
2095 return false;
2096 }
2097 nev=nevt;
2098 }
2099 }
2100 //=========================================================
2101
2102 if( !strcmp(key->GetName(),"NeutronD" ) ){
2103 ND__ok = true;
2104 if(ND){
2105 Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
2106 if( nev && nevt!=nev){
2107 cout << "File: "<< f->GetName() <<" discarded ---- NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
2108 return false;
2109 }
2110 nev=nevt;
2111 }
2112 }
2113 //=========================================================
2114 if( !strcmp(key->GetName(),"Anticounter") ){
2115 AC__ok = true;
2116 if(AC){
2117 Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
2118 if( nev && nevt!=nev){
2119 cout << "File: "<< f->GetName() <<" discarded ---- Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
2120 return false;
2121 }
2122 nev=nevt;
2123 }
2124 }
2125 //=========================================================
2126 if( !strcmp(key->GetName(),"OrbitalInfo") ){
2127 ORB__ok = true;
2128 if(ORB){
2129 Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
2130 if( nev && nevt!=nev){
2131 cout << "File: "<< f->GetName() <<" discarded ---- OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
2132 return false;
2133 }
2134 nev=nevt;
2135 }
2136 }
2137 //=========================================================
2138 if( !strcmp(key->GetName(),"Tracker" ) ){
2139 TTree *T = (TTree*)f->Get("Tracker");
2140 if(TRK1||TRK2||TRKh){
2141 Int_t nevt = T->GetEntries();
2142 if( nev && nevt!=nev){
2143 cout << "File: "<< f->GetName() <<" discarded ---- Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
2144 return false;
2145 }
2146 nev=nevt;
2147 }
2148 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2149 TString name = T->GetListOfBranches()->At(i)->GetName();
2150 if( !name.CompareTo("TrkLevel1") )TRK1__ok=true;
2151 if( !name.CompareTo("TrkLevel2") )TRK2__ok=true;
2152 if( !name.CompareTo("TrkHough") )TRKh__ok=true;
2153 };
2154 T->Delete();
2155 };
2156 //=========================================================
2157 if( !strcmp(key->GetName(),"Calorimeter" ) ){
2158 TTree *T = (TTree*)f->Get("Calorimeter");
2159 if(CAL1||CAL2){
2160 Int_t nevt = T->GetEntries();
2161 if( nev && nevt!=nev){
2162 cout << "File: "<< f->GetName() <<" discarded ---- Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
2163 return false;
2164 }
2165 nev=nevt;
2166 }
2167 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2168 TString name = T->GetListOfBranches()->At(i)->GetName();
2169 if( !name.CompareTo("CaloLevel1") )CAL1__ok=true;
2170 if( !name.CompareTo("CaloLevel2") )CAL2__ok=true;
2171 };
2172 T->Delete();
2173 };
2174
2175 };
2176
2177 if( SELLI==-1 )SELLI = (Int_t)SELLI__ok;
2178 if( SELLI==0 && SELLI__ok ){
2179 cout << "File: "<< f->GetName() <<" discarded ---- found SelectionList (it is not a full-event file)" << endl;
2180 return false;
2181 }
2182 if( SELLI==1 && !SELLI__ok ){
2183 cout << "File: "<< f->GetName() <<" discarded ---- SelectionList missing" << endl;
2184 return false;
2185 }
2186
2187 // cout << "SELLI "<<SELLI<<endl;
2188
2189 // cout<< "CheckLevel2File(TString): detector list --> ";
2190 // if(TRK1__ok)cout<<"TRK1 ";
2191 // if(TRK2__ok)cout<<"TRK2 ";
2192 // if(TRKh__ok)cout<<"TRKH ";
2193 // if(CAL1__ok)cout<<"CAL1 ";
2194 // if(CAL2__ok)cout<<"CAL2 ";
2195 // if(TOF__ok)cout<<"TOF ";
2196 // if(TRG__ok)cout<<"TRG ";
2197 // if(AC__ok)cout<<"AC ";
2198 // if(ND__ok)cout<<"ND ";
2199 // if(S4__ok)cout<<"S4 ";
2200 // if(ORB__ok)cout<<"ORB ";
2201 // cout << endl;
2202
2203
2204 if(TRK2 && TRK1__ok)TRK1=1;
2205 // ----------------------------------------------------------------------------
2206 // NOTA
2207 // se c'e` il level1, lo devo necessarimente leggere.
2208 // infatti (non ho capito perche`) i cluster vengono letti e allocati in memoria
2209 // comunque, ma non vengono disallocati da PamLevel2::Clear()
2210 // ----------------------------------------------------------------------------
2211
2212
2213 if(!RUN__ok) {
2214 cout << "File: "<< f->GetName() <<" *WARNING* ---- Missing RunInfo tree"<< endl;
2215 // return false;
2216 };
2217
2218 if(CAL1 && !CAL1__ok){
2219 cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel1 branch"<< endl;
2220 return false;
2221 };
2222 if(CAL2 && !CAL2__ok){
2223 cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel2 branch"<< endl;
2224 return false;
2225 };
2226 if(TRK2 && !TRK2__ok){
2227 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel2 branch"<< endl;
2228 return false;
2229 };
2230 if(TRK1 && !TRK1__ok){
2231 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel1 branch"<< endl;
2232 return false;
2233 };
2234 if(TRKh && !TRKh__ok){
2235 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkHough branch"<< endl;
2236 return false;
2237 };
2238 if(ORB && !ORB__ok){
2239 cout << "File: "<< f->GetName() <<" discarded ---- Missing ORB tree"<< endl;
2240 return false;
2241 };
2242 if(AC && !AC__ok){
2243 cout << "File: "<< f->GetName() <<" discarded ---- Missing AC tree"<< endl;
2244 return false;
2245 };
2246 if(S4 && !S4__ok){
2247 cout << "File: "<< f->GetName() <<" discarded ---- Missing S4 tree"<< endl;
2248 return false;
2249 };
2250 if(TOF && !TOF__ok){
2251 cout << "File: "<< f->GetName() <<" discarded ---- Missing ToF tree"<< endl;
2252 return false;
2253 };
2254
2255 if(ND && !ND__ok){
2256 cout << "File: "<< f->GetName() <<" discarded ---- Missing ND tree"<< endl;
2257 return false;
2258 };
2259 if(TRG && !TRG__ok){
2260 cout << "File: "<< f->GetName() <<" discarded ---- Missing Trigger tree"<< endl;
2261 return false;
2262 };
2263
2264
2265 // lk->Delete();
2266 // delete lk;
2267 f->Close();
2268
2269 // cout<< "CheckLevel2File(TString): detector list --> ";
2270 // if(TRK1)cout<<"TRK1 ";
2271 // if(TRK2)cout<<"TRK2 ";
2272 // if(TRKh)cout<<"TRKH ";
2273 // if(CAL1)cout<<"CAL1 ";
2274 // if(CAL2)cout<<"CAL2 ";
2275 // if(TOF)cout<<"TOF ";
2276 // if(TRG)cout<<"TRG ";
2277 // if(AC)cout<<"AC ";
2278 // if(ND)cout<<"ND ";
2279 // if(S4)cout<<"S4 ";
2280 // if(ORB)cout<<"ORB ";
2281 // cout << endl;
2282
2283 return true;
2284
2285 };
2286
2287
2288 /**
2289 * Create clone-trees
2290 */
2291 void PamLevel2::CreateCloneTrees0( TChain *fChain, TFile *ofile ){
2292
2293 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2294 cout << "Create clones of PAMELA trees "<<endl;
2295
2296 Int_t i=0;
2297 pam_tree_clone[i] = fChain->GetTree()->CloneTree(0);
2298 TString name = pam_tree_clone[i]->GetName();
2299 name.Append("_clone");
2300 // pam_tree_clone[i]->SetName(name.Data());
2301 cout << pam_tree_clone[i]->GetName() <<endl;
2302 i++;
2303
2304 TList *li = fChain->GetListOfFriends();
2305 TIter next(li);
2306 TFriendElement* T_friend=0;
2307 ofile->cd();
2308 while( (T_friend = (TFriendElement*)next()) ){
2309 // cout<<T_friend->IsA()->GetName()<<" "<<T_friend->GetName()<<hex << T_friend->GetTree() << dec<<endl;
2310 // cout<<T_friend->GetTree()->GetName()<< endl;
2311 pam_tree_clone[i] = T_friend->GetTree()->CloneTree(0);
2312 pam_tree_clone[i]->SetAutoSave(1000000);
2313 name = pam_tree_clone[i]->GetName();
2314 name.Append("_clone");
2315 // pam_tree_clone[i]->SetName(name.Data());
2316 cout << pam_tree_clone[i]->GetName() << endl;
2317 i++;
2318 }
2319
2320 delete li;
2321
2322 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2323
2324 }
2325
2326 /**
2327 * Create clone-trees
2328 */
2329 void PamLevel2::CreateCloneTrees(TFile *ofile){
2330
2331 ofile->cd();
2332
2333 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2334 cout << "Create new PAMELA trees "<<endl;
2335
2336
2337 run_tree_clone = new TTree("Run","PAMELA Level2 data from the GL_RUN table ");
2338 run_tree_clone->Branch("RunInfo","GL_RUN",GetPointerTo("RunInfo"));
2339 cout << "Run : branch RunInfo"<<endl;
2340 run_tree_clone->Branch("SoftInfo","SoftInfo",GetPointerTo("SoftInfo"));
2341 cout << "Run : branch SoftInfo"<<endl;
2342
2343
2344
2345 sel_tree_clone = new TTree("SelectionList","List of selected events ");
2346 sel_tree_clone->Branch("RunEntry",&irun,"runentry/L");
2347 sel_tree_clone->Branch("EventEntry",&irunentry,"eventry/L");
2348
2349
2350 Int_t i=0;
2351 if(TRK1||TRK2||TRKh){
2352 pam_tree_clone[i] = new TTree("Tracker","PAMELA tracker level2 data ");
2353 if(TRK1) {
2354 pam_tree_clone[i]->Branch("TrkLevel1","TrkLevel1", GetPointerTo("TrkLevel1"));
2355 pam_tree_clone[i]->BranchRef();
2356 cout << "Tracker : branch TrkLevel1"<<endl;
2357 cout << "CreateCloneTrees " << GetTrkLevel1()<<endl;
2358 };
2359 if(TRK2) {
2360 pam_tree_clone[i]->Branch("TrkLevel2", "TrkLevel2",GetPointerTo("TrkLevel2"));
2361 cout << "Tracker : branch TrkLevel2"<<endl;
2362 };
2363 if(TRKh) {
2364 pam_tree_clone[i]->Branch("TrkHough","TrkHough", GetPointerTo("TrkHough"));
2365 cout << "Tracker : branch TrkHough"<<endl;
2366 };
2367 i++;
2368 }
2369
2370 // Calorimeter
2371 if(CAL1||CAL2){
2372 pam_tree_clone[i] = new TTree("Calorimeter","PAMELA calorimeter level2 data ");
2373 if(CAL1) {
2374 pam_tree_clone[i]->Branch("CaloLevel1", "CaloLevel1", GetPointerTo("CaloLevel1"));
2375 cout << "Calorimeter : branch CaloLevel1"<<endl;
2376 };
2377 if(CAL2) {
2378 pam_tree_clone[i]->Branch("CaloLevel2","CaloLevel2", GetPointerTo("CaloLevel2"));
2379 cout << "Calorimeter : branch CaloLevel2"<<endl;
2380 };
2381 i++;
2382 }
2383
2384 // ToF
2385 if(TOF) {
2386 pam_tree_clone[i] = new TTree("ToF","PAMELA ToF level2 data ");
2387 pam_tree_clone[i]->Branch("ToFLevel2","ToFLevel2", GetPointerTo("ToFLevel2"));
2388 cout << "ToF : branch ToFLevel2"<<endl;
2389 i++;
2390 };
2391 // Trigger
2392 if(TRG) {
2393 pam_tree_clone[i] = new TTree("Trigger","PAMELA trigger level2 data ");
2394 pam_tree_clone[i]->Branch("TrigLevel2","TrigLevel2", GetPointerTo("TrigLevel2"));
2395 cout << "Trigger : branch TrigLevel2"<<endl;
2396 i++;
2397 };
2398 // S4
2399 if(S4) {
2400 pam_tree_clone[i] = new TTree("S4","PAMELA S4 level2 data ");
2401 pam_tree_clone[i]->Branch("S4Level2","S4Level2", GetPointerTo("S4Level2"));
2402 cout << "S4 : branch S4Level2"<<endl;
2403 i++;
2404 };
2405 // Neutron Detector
2406 if(ND) {
2407 pam_tree_clone[i] = new TTree("NeutronD","PAMELA neutron detector level2 data ");
2408 pam_tree_clone[i]->Branch("NDLevel2","NDLevel2", GetPointerTo("NDLevel2"));
2409 cout << "NeutronD : branch NDLevel2"<<endl;
2410 i++;
2411 };
2412 // Anticounters
2413 if(AC) {
2414 pam_tree_clone[i] = new TTree("Anticounter","PAMELA anticounter detector level2 data ");
2415 pam_tree_clone[i]->Branch("AcLevel2","AcLevel2", GetPointerTo("AcLevel2"));
2416 cout << "Anticounter : branch AcLevel2"<<endl;
2417 i++;
2418 };
2419 // OrbitalInfo
2420 if(ORB) {
2421 pam_tree_clone[i] = new TTree("OrbitalInfo","PAMELA oribital info ");
2422 pam_tree_clone[i]->Branch("OrbitalInfo","OrbitalInfo", GetPointerTo("OrbitalInfo"));
2423 cout << "OrbitalInfo : branch OrbitalInfo"<<endl;
2424 i++;
2425 };
2426 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2427
2428 }
2429
2430
2431 /**
2432 * Fill tree (created with CreateCloneTrees)
2433 *
2434 */
2435 //void PamLevel2::FillNewPamTree(TTree *T){
2436 void PamLevel2::FillCloneTrees(){
2437
2438 // cout << "PamLevel2::FillCloneTrees()" << irunentry << endl;
2439
2440 for(Int_t i=0; i<NCLONES; i++){
2441 if(pam_tree_clone[i])pam_tree_clone[i]->Fill();
2442 }
2443 if(sel_tree_clone)sel_tree_clone->Fill();
2444
2445 }
2446
2447
2448 TTree* PamLevel2::GetCloneTree(TString name){
2449
2450 for(Int_t i=0; i<NCLONES; i++){
2451 if(pam_tree_clone[i]){
2452 TString na = pam_tree_clone[i]->GetName();
2453 if(!name.CompareTo(na))return pam_tree_clone[i];
2454 };
2455 }
2456 if(run_tree_clone){
2457 TString na = run_tree_clone->GetName();
2458 if(!name.CompareTo(na))return run_tree_clone;
2459 }
2460 if(sel_tree_clone){
2461 TString na = sel_tree_clone->GetName();
2462 if(!name.CompareTo(na))return sel_tree_clone;
2463 }
2464 return NULL;
2465
2466 }
2467 void PamLevel2::WriteCloneTrees(){
2468 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2469 cout << "Write clones of PAMELA trees "<<endl;
2470 cout << run_tree_clone->GetName()<<endl;
2471 run_tree_clone->Write();
2472 cout << sel_tree_clone->GetName()<<endl;
2473 sel_tree_clone->Write();
2474 for(Int_t i=0; i<NCLONES; i++){
2475 if(pam_tree_clone[i]){
2476 cout << pam_tree_clone[i]->GetName()<<endl;
2477 pam_tree_clone[i]->Write();
2478 };
2479 }
2480 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2481
2482 }
2483
2484 /**
2485 * Method to get level2-trees entry, the corresponding run entry and (if required) the level0 entry.
2486 */
2487 //Int_t PamLevel2::GetEntry(Int_t iee){
2488 Int_t PamLevel2::GetEntry(Long64_t iee){
2489
2490 if(!pam_tree){
2491 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- level2 trees not loaded"<<endl;
2492 return 0;
2493 }
2494
2495
2496 //
2497 // This is a sort of bug: if you don't have the run tree you don't want to exit here you want to have loaded the event anyway...
2498 //
2499 // if(!run_tree ){
2500 // cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loeaded"<<endl;
2501 // return 0;
2502 // }
2503
2504 Long64_t ii=0;
2505 //-------------------------------
2506 ii = iee;
2507 if( !pam_tree->GetEntry(ii) ){
2508 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- error reading pam tree"<<endl;
2509 return 0;
2510 }
2511 //
2512 // ... that's way I put it here. Notice that nothing change in the code (is backward compatible) since in any case you return with 0.
2513 // in theory one would like to return 1 if run is not loaded but now I don't have the will to add that 2 lines of code and it is not
2514 // a problem if you don't check the return code of getentry.
2515 //
2516 if(!run_tree ){
2517 if ( TRK0 || CAL0 || TOF0 || RUN ) { //forse cosi` va bene per tornare 1?
2518 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loaded"<<endl;
2519 return 0;
2520 } else {
2521 return 1; //cosi` se non c'e` run esce qua...
2522 }
2523 }
2524
2525 //-------------------------------
2526 ii = iee;
2527 // Bool_t UPDATED = UpdateRunInfo(run_tree,ii);
2528 Bool_t UPDATED = UpdateRunInfo(ii);
2529 if(SELLI==0)irunentry = iee-runfirstentry;
2530 if(UPDATED && run_tree_clone)run_tree_clone->Fill();
2531
2532 // cout << "PamLevel2::GetEntry("<<iee<<") "<<irun<<" "<<runfirstentry<<" "<<irunentry<<endl;
2533
2534 if( TRK0 || CAL0 || TOF0 )GetYodaEntry( );
2535
2536 return 1;
2537
2538 }
2539
2540 /**
2541 * Method to retrieve the level0 tree (YODA tree) that contains the current event.
2542 * Given the run ID (...), if needed it query the DB and load the proper file.
2543 * @return Pointer to the tree
2544 */
2545
2546
2547 TTree* PamLevel2::GetYodaTree( ){
2548
2549 // cout << "TTree* PamLevel2::GetYodaTree( )"<<endl;
2550 //===================================
2551 // check if iroot has changed
2552 //===================================
2553 if( irun<0 ){
2554 cout << "PamLevel2::GetYodaTree() -- ERROR -- irun = "<<irun<<endl;
2555 // cout << "In order to use this method you have to first load the RunInfo tree "<<endl;
2556 return NULL;
2557 }
2558 Int_t irootnew = run_obj->ID_ROOT_L0;
2559 // cout << "iroot "<<iroot<<endl;
2560 // cout << "irootnew "<<irootnew<<endl;
2561
2562 //===================================
2563 // load the level0 file
2564 // (if not already loaded)
2565 //===================================
2566 if( iroot != irootnew || !l0_tree){
2567 iroot = irootnew;
2568 //===================================
2569 // open the DB connection
2570 // (if not already opened)
2571 //===================================
2572 if(!dbc || (dbc && !dbc->IsConnected())){
2573 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
2574 cout<<"Connecting to DB"<<endl;
2575 cout<<"HOST "<<host<<endl;
2576 cout<<"USER "<<user<<endl;
2577 cout<<"PSW "<<psw<<endl;
2578 dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
2579 if( !dbc )return NULL;
2580 if( !dbc->IsConnected() )return NULL;
2581 // cout<<"...done"<<endl;
2582 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
2583 }else{
2584 // cout<<"DB already connected"<<endl;
2585 }
2586 GL_ROOT glroot = GL_ROOT();
2587 if( glroot.Query_GL_ROOT(iroot,dbc) ){
2588 cout << "TTree* PamLevel2::GetYodaTree( ) -- ERROR -- level0 file iroot = "<<iroot<< " does not exists"<<endl;
2589 return NULL;
2590 };
2591 TString filename = glroot.PATH + glroot.NAME;
2592 if(l0_file){
2593 l0_file->Close();
2594 l0_file->Delete();
2595 }
2596 cout << "Opening LEVEL0 file: "<< filename << endl;
2597 FileStat_t t;
2598 if( gSystem->GetPathInfo(filename.Data(),t) ){
2599 cout << " PamLevel2::GetYodaTree() -- ERROR opening file "<<endl;
2600 return NULL;
2601 }
2602 l0_file = new TFile(filename);
2603 if( !l0_file )return NULL;
2604 l0_tree = (TTree*)l0_file->Get("Physics");
2605 if(!h0_obj)h0_obj = new EventHeader();
2606 l0_tree->SetBranchAddress("Header" ,&h0_obj);
2607 prevshift = 0;
2608 //---------------------------------------------------
2609 // TRACKER:
2610 if(TRK0){
2611 if(!trk0_obj){
2612 trk0_obj = new TrkLevel0();
2613 trk0_obj->Set();
2614 };
2615 l0_tree->SetBranchAddress("Tracker" ,trk0_obj->GetPointerToTrackerEvent());
2616 }
2617 //---------------------------------------------------
2618 // CALORIMETER:
2619 if(CAL0){
2620 if(!calo0_obj){
2621 calo0_obj = new CaloLevel0();
2622 calo0_obj->Set();
2623 };
2624 l0_tree->SetBranchAddress("Calorimeter" ,calo0_obj->GetPointerToCalorimeterEvent());
2625 // cout << "PamLevel2::GetYodaTree() --- level0 calorimeter not implemented "<<endl;
2626 }
2627 //---------------------------------------------------
2628 // TOF:
2629 if(TOF0){
2630 cout << "PamLevel2::GetYodaTree() --- level0 TOF not implemented "<<endl;
2631 }
2632
2633 };
2634
2635 if(!dbc || (dbc && !dbc->IsConnected())){
2636 cout << " TTree* PamLevel2::GetYodaTree( ) -- no DB connected... hai fatto qualche cazzata "<<endl;
2637 }
2638
2639 if ( TRK0 ){
2640 TrkParams::Load(6);
2641 if( !TrkParams::IsLoaded(6) ){
2642 cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- VK-mask not loaded"<<endl;
2643 };
2644 TrkParams::SetCalib(run_obj,dbc);
2645 TrkParams::LoadCalib( );
2646 if( !TrkParams::CalibIsLoaded() ){
2647 cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- Calibration not loaded"<<endl;
2648 };
2649 }
2650
2651 // cout << l0_tree << endl;
2652
2653 return l0_tree;
2654
2655 }
2656
2657 /**
2658 * Method to retrieve the level0 tree (YODA tree) that contains the current event.
2659 */
2660 Int_t PamLevel2::GetYodaEntry(){
2661
2662 // cout << "Int_t PamLevel2::GetYodaEntry()"<<endl;
2663 if(!GetYodaTree())return 0;
2664
2665 // patch
2666 if( irunentry < 0){
2667 cout << "Int_t PamLevel2::GetYodaEntry() -- ATTENZIONE -- irunentry negativo?!?! "<<(Int_t)irunentry<<endl;
2668 irunentry=0LL;
2669 }
2670 // ---------------------------------
2671 // if file is NOT a preselected file
2672 // ---------------------------------
2673 Long64_t quellagiusta = irunentry + (Long64_t)(run_obj->EV_FROM);
2674
2675
2676 // cout << " irun "<< irun << " irunentry "<< irunentry<<" run_obj->EV_FROM "<<run_obj->EV_FROM <<" quella giusta "<<quellagiusta << endl;
2677 // cout << " iroot "<<iroot<<" run_obj->ID_ROOT_L0 "<<run_obj->ID_ROOT_L0<<endl;
2678 // cout << " time "<< GetOrbitalInfo()->absTime << endl;
2679 // cout << " trk_calib_used "<<run_obj->TRK_CALIB_USED<< endl;
2680
2681 if( !GetOrbitalInfo() ){
2682 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing OrbitalInfo "<<endl;
2683 return 0;
2684 }
2685 if( GetOrbitalInfo()->OBT==0 && GetOrbitalInfo()->pkt_num==0 ){
2686 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- level2 event corrupted ?? "<<endl;
2687 return 0;
2688 }
2689
2690 // ---------------------------------------------------------------------
2691 // ATTENTION!!!
2692 // If data are not pre-processed with cleaner, the level0 tree may contain
2693 // spurious nested physics packets.
2694 // The GL_RUN variables EV_FROM, EV_TO and NEVENTS counts also these entries
2695 // while level2 tree DOES NOT!!
2696 // This means that "quellagiusta" in these cases is not correct.
2697 // In order to retrieve the correct level0 event, I implemented a check
2698 // of the OBT and pkt-number. In case of mismatch, the level0 entry number
2699 // is shift forward until when the packets match.
2700 // ---------------------------------------------------------------------
2701 Int_t answer = 0;
2702 Int_t shift =0;
2703 // printf(" siamo qui %i %i \n",shift,prevshift);
2704 do{
2705 if(shift>0){
2706 cout << " level0 <--> level2 mismatch ( irun "<<irun<<" irunentry "<<irunentry<<" shift "<<shift<<" prevshift "<<prevshift<<" )"<<endl;
2707 }
2708 answer = l0_tree->GetEntry(quellagiusta+(Long64_t)shift+(Long64_t)prevshift);
2709 shift++;
2710 if( !GetEventHeader() ){
2711 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing EventHeader "<<endl;
2712 return 0;
2713 }
2714 // cout << "PKTNUM "<<shift<<" == L2 --- "<< GetOrbitalInfo()->pkt_num << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetCounter()<<endl;
2715 // cout << " L2 --- "<< GetOrbitalInfo()->OBT << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime()<<endl;
2716 // if( (quellagiusta+shift) == l0_tree->GetEntries()+1 )cout << ">>> end of level0 tree <<<"<<endl;
2717 // cout << " GetOrbitalInfo()->OBT "<< GetOrbitalInfo()->OBT << endl;
2718 // cout << " GetEventHeader()->GetPscuHeader()->GetOrbitalTime() "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime() << endl;
2719 // cout << " GetOrbitalInfo()->pkt_num "<< GetOrbitalInfo()->pkt_num << endl;
2720 // cout << " GetEventHeader()->GetPscuHeader()->GetCounter() "<< GetEventHeader()->GetPscuHeader()->GetCounter() << endl;
2721 // printf(" IDRUN %u \n",GetRunInfo()->ID);
2722 //
2723 if ( prevshift != 0 && (quellagiusta+(Long64_t)shift) == GetYodaTree()->GetEntries() ){
2724 prevshift = 0;
2725 shift = -1;
2726 };
2727
2728 }while( ( GetOrbitalInfo()->OBT != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetOrbitalTime()) || GetOrbitalInfo()->pkt_num != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetCounter())) && (quellagiusta+(Long64_t)shift) < GetYodaTree()->GetEntries() );
2729
2730 if ( (quellagiusta+(Long64_t)shift) == GetYodaTree()->GetEntries() ) cout << " Big trouble here, no such event in Level0 data! " <<endl;
2731
2732 // cout << "LA ENTRY GIUSTA E`: "<<quellagiusta<<" (spero...)"<<endl;
2733 // return GetYodaTree()->GetEntry(quellagiusta);
2734 if ( shift > 1 ) prevshift += (shift-1);
2735
2736 return answer;
2737
2738 }
2739
2740

  ViewVC Help
Powered by ViewVC 1.1.23