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

Contents of /PamelaLevel2/src/PamLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.41 - (show annotations) (download)
Thu Aug 2 13:59:47 2007 UTC (17 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.40: +78 -37 lines
track sorting method modified

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

  ViewVC Help
Powered by ViewVC 1.1.23